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Abstract 

A complete formulation is given of an exact kinetic theory for lattice gases. This 
kinetic theory makes possible the calculation of corrections to the usual Boltzmann / 
Chapman-Enskog analysis of lattice gases due to the buildup of correlations. It is shown 
that renormalized transport coefficients can be calculated perturbatively by summing 
terms in an infinite series. A diagrammatic notation for the terms in this series is 
given, in analogy with the diagrammatic expansions of continuum kinetic theory and 
quantum field theory. A closed-form expression for the coefficients associated with 
the vertices of these diagrams is given. This method is applied to several standard 
lattice gases, and the results are shown to correctly predict experimentally observed 
deviations from the Boltzmann analysis. 
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1 Introduction 



Lattice gases are a class of dynamical systems in which particles move on a lattice in discrete 
time steps. In much the same way that the Ising model and other similar lattice spin models 
provide simple examples of equilibrium statistical mechanical systems, lattice gases capture 
many features of nonequilibrium statistical mechanical systems such as fluids. 

Most work on lattice gases has been done in a framework where the Boltzmann assump- 
tion is made, neglecting all correlations between the particles moving on the lattice. In this 
paper we give a complete description of a theory of lattice gases in which the effects of all 
correlations are included. For a wide variety of lattice gases, the effect of these correlations 
is only to correct or renormalize the hydrodynamic transport coefficients, and we show how 
to compute this correction. For some lattice gases, the resulting renormalization can change 
parameters in the theory by substantial amounts, so that in any situation where lattice 
gases are used to make precise quantitative predictions, the effects of correlations should be 
considered, and the magnitude of the resulting correction to the usual analysis should be 
estimated. 



1.1 Background 

In 1986, it was shown Q that lattice gases could be used to simulate two-dimensional Navier- 
Stokes flow. Since then, lattice gases have been developed to describe a wide variety of 
physical systems. Examples are three-dimensional Navier-Stokes flow 0, magnetohydrody- 
namics 0, immiscible fluids with a surface tension interface j|, [5], convection ||, two-phase 
liquid-gas flow @, Burgers' equation |§, and reaction-diffusion equations §|. For a summary 
of recent works on the subject, see the proceedings edited by Doolen [TIJ], by Monaco [fT2"] , 
by Manneville et. al. [f[3| |, by Alves ||14|| , and by Boon and Lebowitz [15 . 



Why should we expect the bulk behavior of particles moving and colliding on a lattice to 
be that of a fluid? In nature, we observe that many different fluids, with drastically differing 
intermolecular force laws, all satisfy the Navier-Stokes equations to a reasonable degree 
of approximation. In spite of all the differences between the intermolecular collisions of, 
say, water and molasses, both types of collisions conserve mass and momentum; ultimately 
the existence of these conserved quantities is what gives rise to fluid-like behavior at the 
macroscopic levelQ. A lattice gas model can be thought of as an attempt to find the simplest 
possible dynamical system with these conservation laws. 

The idea that macroscopic properties of a physical system should be independent of 
the microscopic definition of the system is also a familiar concept in equilibrium statistical 
mechanics and in quantum field theory. In such theories, the effect of looking at the physics 
of the system at larger and larger scales is mathematically described by the renormalization 
group flow of the system [16|]. Generally, as the scale of the physics of interest becomes 



extremely large compared to the scale at which the system is defined, one finds that the 



1 Conservation of mass and momentum are sufficient to get the correct behavior in the incompressible 
regime. To get correct compressible behavior, it is necessary to conserve energy as well. 
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renormalization group flow takes the system towards certain fixed points, which describe 
entire universality classes of theories with identical macroscopic behavior. The emergence of 
similar hydrodynamic equations in a variety of systems with different microscopic dynamics 
is an equivalent phenomenon in nonequilibrium statistical mechanics. 

Of course, there are some properties of fluids in which the details of the microscopic col- 
lisions manifest themselves at the macroscopic level. The most obvious such properties are 
the transport coefficients. While water and molasses both satisfy the Navier-Stokes equa- 
tions, they do so with very different viscosities. Long-time tails in velocity autocorrelation 
functions are also well-known to be sensitive to the microscopic dynamics. 

The central problem in the theoretical analysis of a lattice gas is thus the determination of 
the macroscopic hydrodynamic equations obeyed by the conserved quantities of the system, 
with the transport coefficients expressed as functions of those conserved quantities. This is 
the classical problem of kinetic theory. The methods needed to perform this analysis for con- 
tinuum fluids have been well developed over the past century fT7| . Two key approximations 
are used in the course of such analyses: 

• The Boltzmann molecular chaos approximation (or Stosszahlansatz) neglects correla- 
tions between molecules entering a collision. This makes it possible to obtain a single 
closed equation, called the Boltzmann equation, for the single-particle distribution 
function. 

• The Chapman- Enskog analysis is an asymptotic expansion in Knudsen number^] that 
yields closed hydrodynamic equations for the conserved quantities. 

During the 1960's and 1970's, much work centered on the removal of the first of these two 
approximations for continuum fluids [pl| and for Lorentz gases fl9fl . It was found that the 



buildup of correlations between the molecules of a fluid could seriously alter, or renormalize, 
the transport coefficients predicted by the Boltzmann theory. Expressions for these correc- 
tions were derived in terms of diagrammatic sums. The propagators in these diagrams can 
be thought of as representing correlated quantities, and the vertices as collisions in which 
correlated quantities can interact. Thus, an event in which two particles emerge from a 
collision (and thereby acquire a correlation), move about in a background of uncorrelated 
particles, and later recollide, can be thought of as a one-loop correction to the Boltzmann 
approximation. Further refinements can be obtained by including more intricate diagrams 
- with multiple loops, nested loops, etc. - to account for the interaction of the correlated 
quantities with the background. Standard field-theoretic techniques can then be used to 
approximate these diagrammatic sums. 

As part of the recent flurry of interest in lattice gases, much of classical kinetic theory 
has been extended to lattice gases. In particular, the derivation of Boltzmann equations for 
lattice gases and the application of the Chapman-Enskog theory has been well understood for 
several years now thanks to the works of a number of authors (see, e.g., [20, [2|]). While these 



2 The Knudsen number is the ratio of mean free path to the macroscopic scale length. 
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works have generally capitalized on the similarities between lattice gases and continuum 
fluids, they have also pointed to some very important differences between the two. 

The most important of these differences has to do with the assumptions leading up to 
the derivation of the Boltzmann equation. Though the assumption of molecular chaos is 
made in either case, the analysis of continuum fluids generally proceeds under the additional 
assumption that the fluid is dilute - that is, that collisions involving more than two molecules 
can be ignored. For lattice gases, on the other hand, it is essential that tertiary and higher 
collision events be treated properly. For example, tertiary collisions are essential to the 
success of the FHP-I lattice gas because they break an unphysical (spurious) conserved 
quantity In addition, the inclusion of higher-order collision events has been shown to 
decrease the viscosity of fluid lattice gases, thereby making possible higher Reynolds' number 
simulations [ED|. 



Thus, collision operators for lattice gases routinely include terms for tertiary and higher- 
order collisions. While this is really only a minor nuisance in the Boltzmann and Chapman- 
Enskog analyses, it introduces a significant complication to the exact kinetic theory. Con- 
tinuum kinetic theory needs to consider only those vertices for which two or fewer correlated 
quantities enter and exit. Lattice gas kinetic theory, on the other hand, must additionally 
treat vertices involving three or more entering and/or exiting correlated quantities. The 
number of correlated quantities that can enter and/or exit a vertex is limited only by the 
number of lattice vectors at a site. This gives rise to a much richer diagrammatic series than 
is the case for a continuum fluid. 

Another difference between continuum fluids and lattice gases, which gives the kinetic 
theories for these types of systems very different flavors, is the simple fact that the contin- 
uum fluids are defined on continuous spaces while the lattice gases are defined on discrete 
spaces. The consequence of this difference is that while the corrections due to correlations 
in continuum fluids are expressed in terms of complicated integral expressions, the correc- 
tions in lattice gases are given by combinatorial sums over countable sets of graphs. Thus, 
the problem of calculating exact transport coefficients in a lattice gas becomes essentially a 
combinatoric, rather than an analytic, problem. 

For some lattice gases the dominant corrections to the transport coefficients arise from 
diagrams which extend spatially and temporally over only a short distance on the lattice^]; 
this dramatically simplifies the numerical calculation of renormalized transport coefficients. 
The resulting simplicity in the formulae describing diagrammatic corrections to transport 
coefficients in lattice gases make these systems extremely attractive both as models with 
which to study complicated properties of related physical systems in the continuum, and as 
pedagogical tools with which to describe the essential features of a complete kinetic theory 
without the complications of many continuum systems. 



i.e., a distance that goes to zero in the scaling limit 
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1.2 Results Presented in This Paper 

In this paper, we give for the first time a complete formulation of an exact kinetic theory of 
lattice gases. We prove that for a wide class of lattice gases, the effects of correlations can be 
completely described by a renormalization of the transport coefficients. Furthermore, we give 
explicit formulae for the renormalized transport coefficients in terms of a sum over an infinite 
set of diagrams. We describe the renormalization calculation in detail for several specific 
lattice gases, and discuss a variety of methods for approximating the complete diagrammatic 
sum. 

We begin with a review of the Boltzmann and Chapman-Enskog theories. Though, as 
has been noted, these methods have been successfully extended to lattice gases for several 
years now, our presentation contains several novel features and results: 

• Our final results for the hydrodynamic equations and the transport coefficients have, 
to our knowledge, never before been written down in this generality. 

• Our results are presented in a mathematical format that is interesting in its own right. 
We introduce a family of metrics on the space of the Boltzmann distribution function, 
and show that the transport coefficients can be expressed quite naturally in terms of 
these metrics and their associated Christoffel symbols. 

• For the special case of lattice fluids on a regular lattice, we derive an expression for 
the viscosity that agrees with that of Henon, and yields a new and compact form for 
the quantity that he calls A in his paper on the subject pifl . 

Next we present the exact kinetic theory for lattice gases. The main features of this 
theory are as follows: 

• We give a simple expression for the renormalized transport coefficients in terms of 
a sum over an infinite set of terms, each associated with a diagram describing the 
propagation and interaction of a set of correlations between lattice sites. 

• We show that the correction associated with each diagram can be decomposed into a 
product of elementary factors associated with the vertices in the diagram; these factors 
describe the interactions of correlated quantities at a single lattice site. Furthermore, 
we give a closed-form expression for these vertex coefficients that is valid for any lattice 
gas obeying semi-detailed balance at lowest order in the Chapman-Enskog expansion 
parameter. 

• We show that when a lattice gas violates semi-detailed balance at higher order in the 
expansion parameter, the effects of correlations can give rise to renormalized source 
terms in the hydrodynamic equations in addition to the usual renormalization of trans- 
port coefficients. We describe these renormalized source terms in terms of the same 
types of diagrams used to describe the renormalized transport coefficients. 
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We formulate the theory in terms of connected correlation functions. A similar expan- 



sion in terms of products of fluctuations was previously used in |22], |23[ and in |24], p5 
where the ring kinetic theory was derived for lattice gases, and for lattice Lorentz gases, 
respectively. The use of the connected correlation functions simplifies the form of the 
complete diagrammatic expansion for the kinetic theory significantly. 

We describe a variety of simplifications and approximations which can be made to 
simplify the numerical computation of the renormalized transport coefficients. The 
simplest of these approximations is the ring approximation, in which correlations be- 
tween more than 2 particles are neglected. Thus, we find that the results of []22| are 
described in our theory by the restriction of the diagrammatic sum to the simple set 
of ring diagrams. 

Previous work on the kinetic ring theory of lattice gases has used the Green-Kubo 
formalism to obtain the series for the transport coefficients. In this paper, we use the 
Chapman-Enskog theory instead. In this way, we get the Boltzmann approximation at 
zeroth order. All higher-order terms in our series are thus corrections to the Boltzmann 
approximation. In the Green-Kubo theory, by contrast, it is necessary to sum an 
infinite number of terms just to get the Boltzmann approximation. 

We work out the diagrammatic expansion in detail for four model lattice gases, and 
we compare some of these results to lattice gas simulations. 

For one of these model problems - a lattice gas for Burgers' equation || - we show that 
the Boltzmann approximation is exact. This result has been previously obtained by 
other means p6| . We show that the diagrammatic formalism introduced here reduces 



the proof to essentially a triviality. 
1.3 Layout of This Paper 

Section ^| begins with a review of lattice gas theory, and establishes notation that is used 
throughout the rest of the paper. Section |3] then reviews the Boltzmann and Chapman- 
Enskog analyses, and introduces a novel mathematical formalism with which the results for 
the transport coefficients can be expressed in a simple form. Section |] contains examples in 
which we work out the details of this formalism for four model lattice gases. 

Section |5| formulates the complete lattice gas kinetic theory, proves that the exact trans- 
port coefficients are expressible as the sum of a diagrammatic series, and gives a closed-form 
expression for the vertex coefficients. Section ^| then works out these vertex coefficients for 
the same four model lattice gases described in Section |], and presents the results for the 
transport coefficients of these lattice gases as formal diagrammatic series. For the Burgers' 
equation lattice gas we show that the diagrammatic series of corrections to the Boltzmann 
approximation vanishes, and thereby recover the result of . 

Since, with the exception of the Burgers' equation example, these diagrammatic series are 
difficult to sum without resorting to numerical methods, Section [5] includes a discussion of 
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approximation schemes for these diagrammatic sums, including summations for short times, 
summations over one- loop diagrams (ring approximation), and summations over subsets of 
diagrams which correspond to truncations of the BBGKY hierarchy. In Section [7], we apply 
several of these approximation schemes to one of the model problems (the 1D3P lattice gas), 
and we present comparisons to experiment, which show that our theory correctly predicts 
observed deviations from the predictions of the Boltzmann theory to within measurable 
accuracy. In the course of doing this, we expand on the discussion of approximations in 
Section [|, further developing techniques valid for general lattice gases. 

At the end of the paper, a glossary of notation is included due to the large number of 
symbols introduced in the text. 



7 



2 Review of Lattice Gas Theory 



2.1 Definitions and Notation 

A lattice gas is generally described by a state space and a time-development rule. The state 
space is defined by associating n bits with each point on a lattice L. (Bits are variables 
taking values in {0, 1}.) We define the set of bits at a general lattice site to be B, so that 
\B\ = n. We denote the total number of bits on the lattice by iV = n\L\. For each value of 
the discrete time parameter t, we write the values of the bits as n*(x, t), where xel and 
i G B. As noted in Section [j], these bits can be thought of as a set of occupation numbers 
for individual particle states. 

We define the set S of possible states of the bits at a general lattice point at a fixed value 
of t to be 

S = {s : s C B}, 

where a state s is associated with the set of bits taking the value 1. We write the value of 
bit i in state s as 




1 if i G s 
otherwise 



Note that \S\ = 2 n . 

We shall often want to refer to the N = n\L\ bits of the lattice in a uniform fashion, so 
we introduce an enumeration of these N bits, given by a 1-1 correspondence between the 
sets B = {1,2, . . . ,N} and B x L = {(z,x) : i G 5,x G L}. In this notation, a single bit 
of the lattice gas is written as n a , a G B. To relate this notation to the more explicit (i,x) 
notation, we express the above 1-1 correspondence by writing a and (i, x) as functions of 
one another, so that 

n a (t) = n i(a \x(a),t) 

and 

n 4 (x,t) =n a{t ' x) (t). 
We shall use both notations interchangeably throughout this paper. 

2.2 Microscopic Dynamical Equation 

For each value of i G {1, . . . , n} there is a lattice vector c l such that x + c l G L for every 
x G L. The evolution of a lattice gas for one timestep can be divided into two substeps: (1) 
A collision substep in which the n bits at each site may alter their values to model a local 
interaction between the particles that they represent, and (2) an advection substep in which 
the new value of bit i at site x moves to that of bit i at site x + c\ 

We wish to express this evolution in the form of an equation for the microscopic dynamics 
of the lattice gas; that is, we desire an equation for n l (x, t + At) in terms of n l (x, t), where 
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At denotes the timestep. Suppose, for a moment, that the particles simply advected without 
colliding. Then the dynamics would be described by 

n*(x + c\ t + At) = n i (x,t). 

The addition of collisions introduces a collision operator on the right-hand side of the above 
equation. That is, we have 

n i (* + c\t + At) =n i (x,t) + c i (n*(x,t)), (2.1) 

where the collision operator, c l , describes the change in bit i due to collisions. Note that 
the form f(z*) is used to indicate that a function / depends on a quantity z for all possible 
values of the index that has been replaced by the asterisk. 



2.3 The Collision Operator 

The collision process at a fixed lattice site and timestep can be fully specified by a 2 n by 2 n 
Boolean transition matrix, a, whose element a(s — > s') is unity if and only if the particles in 
state s collide to yield particles in state s'. Since each incoming state gives rise to exactly 
one outgoing state, 

J2a(s^s') = l. (2.2) 

s' 

Stochastic lattice gases deserve special mention at this point. They also have a(s — > s') 
equal to unity for exactly one value of s' for each value of s, but that value of s' may differ 
from site to site, and from timestep to timestep for a fixed value of s. For example, we might 
have a(s — ► s±) = r, a(s — > S2) = 1 — r, and a(s — > s') = for all other s', where r is a 
random bit that is sampled at each site at each timestep with some specified mean, R — (r) . 
This would effectively mean that the outcome of a collision for incoming state s is state s\ 
with probability R, and state S2 with probability 1 — R. Note that Eq. Q2.2p is still always 
satisfied. 

If we define the Kronecker delta function of two bits , 



8{x,y) 

then the product 




(2.3) 



3=1 

is unity at a given site if and only if that site is in state s. Since the collision operator c l is 
nothing more than the total change in bit % due to collision, it can be expressed in terms of 
the transition matrix as follows: 

n 

c\ n *) = -> s')(s H - s*) J] 8(n j , s j ). (2.4) 

s,s' j=l 
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2.4 Conserved Quantities 

Another distinguishing feature of lattice gases is the presence of some number of additive 
conserved quantities that are linear in the bit values. For example, in some lattice gases, 
the total number of particles is a conserved quantity; it is clearly conserved by the advection 
phase of the timestep, and we can choose collision rules that conserve particle number as 
well. 

Let us assume that we have a lattice gas with some number, n c , of conserved quantities. 
We assume that all conserved quantities are linear in the configuration bit values, so we can 
write the value of the /ith conserved quantity at site x and at time t in the formP] 

g^(x, t) = gf n*(x, t) for fi = l,..., n c . (2.5) 

where the coefficients gf satisfy 

= gf c \n*) for n = 1, . . . , n c . (2.6) 

It is important to note that the gf are constant coefficients, independent of spatial position. 
From Eq. ( |2.6|) , it follows that 

£g^(x,t + At) = £g<V(x + c\t + At) 

X X 

= ]Tgf [n'(x,t)+cV)] =£g"(x,t), 

X X 

so that the global sum of any conserved quantity is constant in time. 

Note that the conserved quantities naturally partition the set S into equivalence classes. 
Two states belong to the same equivalence class if they have the same values of all the 
conserved quantities (that is, s ~ s' iff gfs /l = gfs* for \i = 1, . . . ,n c ). Collisions must map 
states into other states of the same equivalence class. This means that the transition matrix, 
a(s — > s') is block diagonal, in that a(s — > s') = if s ^ s'. More succinctly, 



a(s 



Of - tf**) = 0. (2.7) 



Note that Eq. ( ^.6[ ) follows immediately from Eqs. ( |2.4j ) and ( J2 . 7| ) . 

We observe in passing that some lattice gases also posess spurious global conserved quan- 
tities. For example, it is easily seen that lattice gases of single-speed particles on a Cartesian 
grid conserve all quantities separately on both checkerboard sublattices. Such spurious global 
conserved quantities have no analog for continuum fluids, and need to be considered carefully 
when using lattice gases to model hydrodynamic phenomena. (See, e.g., f29||.) 



4 Throughout this paper we adopt the modified Einstein summation convention that if the same index 
appears in at least one contravariant position and at least one covariant position in every term where it 
appears at all, then it is to be summed over its entire range of values. Sometimes summations are indicated 
explicitly, particularly when the range of summation is not clear from the context. 
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2.5 The Ensemble Average 



We now consider some statistical aspects of lattice gas theory. Let us suppose that we 
have prepared an ensemble of lattice gas simulations, on grids of the same size, with initial 
conditions that are identical at the hydrodynamic level but may differ (are sampled from 
some distribution of initial conditions) at the kinetic level. We may then take averages across 
this ensemble. We shall denote these ensemble averages by angular brackets (). We denote 
the ensemble average of the quantity n l (x, t) by iV l (x, t); that is 

AP(x,t) = (r/(x,t)> 

Note that while the n*'s are binary, the N l, s take their values in the set of real numbers 
between zero and one. 

Next, the ensemble average of a(s — > s') is defined to be 

A(s -> s') = (a{s -> s')), 
so that the ensemble average of Eq. Q2.2|) is 

£A(* -✓) = !■ (2-8) 

s' 

For deterministic lattice gases, A(s — > s') = a(s — > s'), whereas for stochastic lattice gases, 
the elements of A(s — * s') are generally real numbers between and 1. 

Similarly, we can consider ensemble-averaged values of the conserved quantities 

Q"(x, t) = (g"(x, £)> = q? A^x, t) for = 1, ... , n c . (2.9) 

Note that we now have three possible levels of description of the lattice gas system. At 
the finest level, the specification of the n l (x, t) constitute a complete microscopic description 
of the system. Their time evolution can generally be obtained only by an actual simulation 
of the lattice gas. Often, however, precise knowledge of each and every bit of the system is 
more information than one really desires. A coarser description, such as a closed set of kinetic 
equations for the ensemble averaged N l (pc,t) is often a more appropriate description of the 
system. Even this level of description, however, is redundant for many purposes. Therefore, 
at the coarsest level, one might seek a closed set of hydrodynamic equations for the The 
remainder of this paper will be concerned with deriving these two reduced descriptions of 
the system. 

2.6 The Boltzmann Equation 

Toward the goal of obtaining a closed set of equations for the iV*, we take the ensemble 
average of Eq. (|2.1j). We are immediately thwarted by the fact that the collision operator, 
given by Eq. Q2.4|) , is generally a nonlinear function of the n % . As is well known, the average of 
a nonlinear function is not in general expressible as a function of the averaged quantitites. It 
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also depends on the correlations between the quantities - in this case between the incoming 
bits, ri 1 . 

Thus, the simplest approximation that we can make to close the system of equations 
for the N l is to assume that the incoming bits, n % , are uncorrelated. This is the discrete 
version of the famous Boltzmann molecular chaos assumption. From this assumption, it 
would follow that 

{din*)) ^C i ((n*)) = C i (N*), 



where, using Eq. ( fOD and Eq. (fO|), C l (N*) is given by 

n 

C\N*) = J2 A (s^ s')(s H - s*) J] (l - s J - W + 2s j N j ) 

s,s' j=l 

= J2 A (s^s')(s' i -s i )f[ (N j ) s3 (l-N j y~ S \ (2.10) 

s,s' j=l 

In this way, we get the lattice Boltzmann equation, 

N i (x + c\t + At) =iV i (x,t) + C i (iV*(x,t)). (2.11) 

Physically speaking, the assumption of molecular chaos supposes that the advection 
substep effectively decorrelates the different bits at each siteQ. That is, it supposes that 
colliding particles have never had any prior effect on each other. This assumption is virtually 
never strictly correct for a system of particles moving on a discrete lattice in a finite number 
of dimensions. By standard combinatorial arguments, the reencounter probability for two 
particles executing a random walk on a lattice is unity in one and two dimensions, is less 
than unity in three or more dimensions, and falls to zero as the number of dimensions goes 
to infinity. One might thus expect that the molecular chaos assumption becomes more valid 
as the number of spatial dimensions increases. Indeed, this is the case, and the molecular 
chaos assumption can be thought of as a sort of mean-field theory. In addition, in some 
circumstances, it is possible for particles to set up coherent structures that persist for long 
times. Such structures, by their very nature, invalidate the molecular chaos assumption in 
a rather dramatic way. 

The remainder of this paper is devoted to deriving the desired closed set of hydrodynamic 
equations for the Q 11 . We shall go about this task in two stages. First, in Section [5], we shall 
show how they can be derived under the molecular chaos assumption. In Section |5|, however, 
we shall abandon this assumption. We shall find that this has the effect of correcting, or 
renormalizing, the transport coefficients in the resulting hydrodynamic equations. For a 
large class of lattice gases of interest - those satisfying a condition known as semi- detailed 
balance at lowest order - we shall show that it is possible to write an exact expression for 
this correction diagrammatic series. 



5 For stochastic lattice gases, such decorrelation is enhanced by the injection of stochasticity at each site 
at each time step. 
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3 The Chapman- Enskog Analysis 



3.1 Asymptotic Ordering 



We shall now outline a perturbative analysis of the Boltzmann equation, Eq. ( [2.11Q . With 



this analysis, we study the hydrodynamics of systems which deviate slightly from local 
equilibrium conditions. To do this, we must first establish the asymptotic regime that we 
are trying to study In this paper we use what is sometimes called diffusion ordering or 
Navier-Stokes ordering. This ordering can be obtained formally by letting c — > ec and 
At — ► e 2 At in the dynamical equations, where e is an expansion parameter. Thus, we are 
taking At ~ c 2 , as is appropriate for diffusive or viscous processes. 

Because the N l are real numbers (as opposed to the n l which are bits), we are free to 
approximate them by smooth functions that happen to coincide with them in value on the 
lattice points. We can then Taylor expand Eq. ( |2.11| ), retaining terms to order e 2 . We get 

3 W e 2 
e 2 At-^- + ec* ■ VN { + -cV : VViV* = C l {N*), +eC{{N*) + e 2 C l 2 (N*), 

where the double-dot notation (:) denotes two inner products^. Note that we have ordered 
the ensemble-averaged collision operator in the expansion parameter e. In what follows, we 
shall assume that Co and C\ respect the conservation laws exactly, but that C 2 does not 
necessarily do so. This will allow us to consider lattice gases whose conservation laws are 
only approximate. 

The above equation can be written in the more suggestive form, 

f)N i ( v i \ ( rV' \ 1 f f 1 

' V + £V ■ (a^J + £2VV ^ {m N ') = Ai^ + Ai C ^> + A^(A"). (3.D 

By contracting this with the constant gf , we get the n c conservation equations 




^Ci(N*), (3.2) 



where /x = l,...,n c . We can now clearly identify the quantity in brackets as the flux 
corresponding to the conserved density and the right-hand side as a source/sink term. 

In what follows, we shall expand the N l in a perturbation series in powers of e about an 
equilibrium state, 

JV* = N* + eN[ + e 2 N l 2 + ■■■ (3.3) 

Here Nq is a local thermodynamic equilibrium. In the next section, we shall characterize 
these equilibria. Then we shall derive hydrodynamic equations for the system by considering 
its near-equilibrium behavior. 

6 This notation will be used only when no ambiguity can arise from it, as is the case when at least one of 
the two dyads involved is symmetric 
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3.2 Semi-Detailed Balance and Equilibria 



An equilibrium distribution for a given lattice gas is a distribution on the state space of the 
system that is invariant under the full dynamics. A Boltzmann equilibrium is a set of values 
for the mean occupation numbers, N a , which is invariant under the Boltzmann equation, 
Eq. ( |2.11D . A Boltzmann equilibrium can be associated with a distribution on the set of states 
by independently sampling each bit n a with probability N a . If the Boltzmann equilibrium is 
spatially uniform, it can be specified by a set of values for the n mean occupation numbers, 
N l ] the associated distribution on the set of states is given by independently sampling each 
bit n a with probability N^ a \ 

A Boltzmann equilibrium is defined to be stable when it is also an equilibrium in the 
more general sense; that is, when the full dynamics of the system, Eq. ( j2.1|) , do not generate 
correlations between the n a . 

We wish to study the dynamics of the system in the vicinity of a local equilibrium. Thus, 
we demand that the lowest-order terms in the mean occupation numbers, Njf, correspond to 
a Boltzmann equilibrium with respect to the lowest-order part of the Boltzmann equation, 
Eq. ( |3.1|) . This condition is given by 

Cg(JVJ(x,t))=0. (3.4) 

We will also demand that the zero-order equilibrium Nq be stable in the sense that the 
zero-order collision operator does not generate correlations between the bits of the system. 
However, we shall allow the zero-order equilibrium parameters iVg(x, t) to have spatial de- 
pendence. In this paper, we shall henceforth restrict our attention to lattice gases for which 
equilibria with these desired properties exist. Fortunately, the restriction that this places on 
the lattice gases that we may consider can be stated as a simple sufficient condition on the 
lowest order part of the transition matrix [Q. To state this condition, we first introduce a 
definition: 

Definition 1 A lattice gas is said to obey detailed balance if its transition matrix satisfies 

A{s -> s') = A{s' -> s), (3.5) 
and it is said to obey semi-detailed balance if its transition matrix satisfies 

J2A(s^s') = l. (3.6) 

s 

Note that semi-detailed balance is a weaker condition than detailed balance, because 
Eqs. (|3~5| ) and ( |2.8| ) together imply Eq. (|3.6| ). Semi-detailed balance (coupled with probability 
conservation) requires that the rows and columns of the transition matrix all sum to unity. 
Detailed balance additionally requires that it be symmetric. 

For a deterministic lattice gas, detailed balance requires that if the transition matrix 
takes state s to state s', then it must also take state s' to state s. In this case, semi-detailed 
balance is the weaker condition that the final states are a permutation of the initial states. 

We now quote [Q] and prove a theorem on the existence of stable Boltzmann equilibria. 
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Theorem 1 Stable and spatially uniform Boltzmann equilibria exist for any lattice gas obey- 
ing semi- detailed balance. These equilibria are described by the Fermi-Dirac distribution, 



No = ^ (3-7) 



l + exp(-E^a a^q? 



where the are n c arbitrary multipliers. 



Proof: At first, Eq. ( |3.4| ) appears to impose n conditions on the n unknowns, Nq, but the n c 
restrictions imposed by Eq. (|2.6|) mean that only n — n c of these conditions are independent, 
and therefore that we ought to expect an n c parameter family of equilibria. These will be 
the constants, a M in Eq. ( |3.7| ). 



Taking the Fermi-Dirac distribution ([3.71) to define an independent distribution on the 



bits of the system, the probability at any site of a fixed state s is given by 

Po(s) = n (my J (i - mf~ s \ ( 3 .8) 

i=i 

In order to show that this distribution on states defines a stable equilibrium, we must show 
that the distribution is unchanged by the collision operator. From the block diagonal prop- 
erty ( |2.7| ) of the transition matrix, and the property of semi-detailed balance, it is clear that 
it will suffice to prove that the probability (|3.8|) is dependent only on the equivalence class 
of s, and therefore on the quantities gf s l for /i = 1, . . . , n c . We have 

„ , , ex P (- ££=i a„ Ei QiO- ~ s 
p o(s) ~ 



Ui [l + exp(- E;=i « M gf 



This expression indeed only depends upon the equivalence class of the state s, and therefore 
we have proven that the distribution on states defined by Eq. ( p.7|) is a stable equilibrium. 

□ 

The equilibrium distribution given by Eq. (|3.7|) is an n c -parameter family of solutions for 
the Nq. By taking the parameters a M to be spatially dependent, we can construct a family 
of Boltzmann equilibria which are not equilibria in the more general sense, but which are 
still stable under the collision operator. These are precisely the type of spatially varying 
equilibria which we desire for the lowest order means Nq of our lattice gas. Thus, in what 
follows we restrict our attention to lattice gases that obey semi-detailed balance at lowest 
order. 

Summarizing the constraints on the collision operator at each order, we have the following: 
C must respect the conservation laws and obey semi-detailed balance; G\ must respect 
the conservation laws, but may violate semi-detailed balance; C*2 can violate either the 
conservation laws or semi-detailed balance. 

Finally, we note that since the parameters, a^, are arbitrary multipliers, any set of 
n c independent functions of them would also suffice to parametrize the equilibrium. In 
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particular, a natural and logical choice of parameters are the hydrodynamic densities, Q^. 
These can be related to the a M by their definition (see Eq. ( |2.9| )), 



Q M = q?N* (a*) for l i= l,...,n c . (3.9) 

Thus, the equilibrium distribution, Nq can be parametrized solely by the equilibrium values 
of the n c conserved densities. 

3.3 The Fermi Metric 

In what follows, we shall need the first two derivatives of the Nq with respect to the so 
we compute them here by the chain rule. First, by differentiating Eq. ( |3.7| ) with respect to 
a u , we obtain 

£-£ = JV$ (l - N$) q\ (3.10) 

(where there is no summation over % because it appears only once on the left). Next, by 
differentiating Eq. ( |3.9| ) with respect to Q u , we obtain 

v Pi <H d Q v PC d Q"' 

where we have defined the symmetric rank- two tensor, 

gX = Nl(l-Ni)tf& (3.11) 
We denote the inverse of this matrix by g^ v so that 

n c 

Since g^ u is a symmetric second-rank tensor, we can identify it as a metric on the space of 
hydrodynamic variables, Q u . We call it the Fermi metric. In terms of the Fermi metric, we 
have 

~dQ v = 9iv ' 

Finally, we can write 



(7 A 7, | 

v=i 



Q^ = EM 1 ~ N i)^ (3 - 12) 



(no sum on i). 

In similar fashion, we compute the second derivative 

°- = JV*(1 - N*)(l - 2Ni)qUl9^ + 2^(1 - N^g^T\ u , (3.13) 



8Q^dQ h 
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where we have defined the Fermi connection, 



2 y \dQ» dQ» dQt 



<Ni (l - W) (l - 2Ni) q]qUl 



(3.14) 



Next, we introduce a characteristic lattice spacing c to define the dimensionless lattice 
vectors, 



i c 

c 



and consider the completely symmetric outer product of k of these vectors, <8> fe e\ It will be 
useful to include these outer products in the above sums. Thus, we define the generalized 
Fermi metric, 



S{k)* = N' (l - N$) qfql 
and the generalized Fermi connection , 

1 



r(*)' 



/J.U 



g^Nt (l - m) (l - -_>.V() q)q]q] 



(3.15) 



(3.16) 



Note that g(0) M ^ = and r(0) T? = Once again, to raise and lower the indices of 

these objects, we use the Fermi metric, g uti , as a metric tensor; thus, e.g., g(2)^ = g(2)^p Atl/ . 
In passing, we note that since the e J are obviously independent of the Q*, we have 



r(ky 



\9 



'dg(k) 



dg(k) 



for all k. That is, we have defined a set of metrics and their associated connections. For 
each k, the members of this set comprise a completely symmetric tensor of rank k. As we 
shall show, this structure is very useful for the problem at hand. 
From Eqs. ( gig) and (|3~T3| ), we get 







and 



d 2 



OQ^dQ 1 



qlK 



q% iv o 



Note that when k = these reduce to the identities dQ^/dQf 1 = ^ and d 2 Q^/dQ^dQ v = 0, 
respectively. 
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3.4 Zero-Order Conservation Equations 

We can now examine the conservation equation, Eq. (|3.2|) , at 0(1). We have 



V- [q?hK =° for/x = l,...,n c 



(3.17) 



Using Eq. ( |3.12 ), this can be written solely in terms of the conserved densities and their 
gradients. We find 

sEI^gayVW^O for // = !,..., n c . 



3.5 The Linearized Boltzmann Equation 

We return to Eqs. Q and Q. At 0(e) we find 

1 



V • — iVn - — 

[At At 



J\N{ + C\{Nq) 



(3.18) 



where we have defined the Jacobian matrix of the lowest-order collision operator at equilib- 
rium, 



3 ~ dW 



N=N 



By differentiating Eq. ( |2.10| ), we can write this directly in terms of the lowest-order transition 
matrixF], A Q (s — > s 1 ), 



J) = EM" - s')(s H - s%2s k - 1) n «f (l - N ( 



l-s k 



(3.19) 



Note that the gf s comprise the components of n c null left eigenvectors of J'-, since 

d 



^ p. 



0. 



N=N 



In what follows, we denote the eigenvalues of J by A M . The corresponding right (left) 
eigenvectors are indexed as contravariant (covariant) vectors, and enumerated with a sub- 
script (superscript). Thus, 

J)qi = \% (3.20) 

and 

q?J) = A"#. (3.21) 

Note that, for = 1, ...,n c , this coincides with the definition of the gf introduced in 
Eq. (|3). 



7 That is, the transition matrix corresponding to the lowest-order collision operator, Cq. 
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The modes enumerated 1, . . . ,n c , correspond to null eigenvalues of J, and will be called 
hydrodynamic modes; those modes enumerated n c + 1, . . . , n will be called kinetic modes. 
We shall often write H for the set {1, . . . , n c }, and K for the set {n c + 1, . . . , n}. For lattice 
gases of interest, the kinetic eigenvalues satisfy \ u < 0; in this case we define the system to be 
linearly stable. Since the eigenvalues of J set the time scale for the equation, we see that for 
a linearly stable lattice gas, the kinetic modes decay away rapidly, while the hydrodynamic 
modes persist for long times. For the remainder of this paper we restrict attention to lattice 
equilibria which are linearly stable. 

Postmultiplying Eq. (|3.21| ) by q^ premultiplying Eq. (|3.20|) by gf , and subtracting, we 

get 

= (X 13 - y)qVq% 

so that right and left eigenvectors corresponding to different eigenvalues are orthogonal. 
Thus, they may be chosen so that 

= (3-22) 



By including (g) fc e* in Eq. ( ^T22|) , we can define a generalized Kronecker delta in the same 



spirit that we generalized the Fermi metric and connection. Thus, 

W>9?((8)e J 'W (3.23) 



so that <5(0)^ = 5^ u . Also, note that the indices of the generalized Fermi metric and 
connection, as well as those of the generalized Kronecker delta, can now be extended to run 
over the kinetic modes as well as the hydrodynamic modes, by simply using the kinetic left 
eigenvectors in Eqs. ( |3.15| ), ( [3.16|) and (|3.23|) , respectively. 

Finally, we note that it is possible to write an explicit, closed expression for the right 
hydrodynamic eigenvectors. From Eqs. (|3~^1), (|3~7|), and (|3.10|) , it follows immediately that 



= fofoW) = J)q£ = 4 [^(1 - JVg) 
whence we identify the right hydrodynamic eigenvectors, 

ti = <$N*(l-Nt). 

3.6 First-Order Solution 



Consider Eq. (|3.18|) for the A 7 ]. Since J is a singular matrix, we must verify that the equation 
is consistent. The consistency requirement is found by premultiplying the equation by the 
null left eigenvectors. We obtain the requirement 




±tf (J^NI + Cl(N*)) 
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for \i G H. Note that this consistency requirement is precisely the zeroth-order conservation 
equation, Eq. ( |3.17| ). It follows that V ■ (c'JVJ) — C[(Nq) has no components in the null space 
of J, since we have assumed linear stability. Thanks to the completeness and orthonormality 
of the eigenvectors, this expression can be written as 



where 



if 



V • (tf'JVg) - C{(N*)] for ueK. 



The solution for N[ can then be written down immediately 

where the 9 U are arbitrary. We fix the solution by assuming, without loss of generality, that 
Q v = 0. That is, we assume that the first order contribution does not affect the definitions 
of the conserved densities. The final result for the N[ is then 



QlQj 
X v 



C((N* ) 



so that we have 



N o + * E 4# [V • (c^) - C((N*)} + G(e 2 ). 



X" 



(3.24) 



(3.25) 



3.7 First-Order Conservation Equations 

We now write the conservation equation, Eq. ( ^2|) , retaining terms to 0(e). We get 

'(tf&W) + V • (?f fgiV*)] = &tfC*(N5) for li e H. 



at ^ v 



Substituting Eq. (|3.24 ) for the N[, we get, after some manipulation, the following closed set 
of equations for the conserved densities: 



^ + V • A» 



E^h V ■ (X>^ ■ VQ*J + S» for fx EH, 

where the advection, diffusion and source coefficients are given by 



WW) - ^7 E 



At 



E 



i 



;(2/ 



(3.26) 

(3.27) 
(3.28) 
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S»(Q*) = ^(Q*), (3-29) 

respectively, and where 

CH(Q*) = q^Ci(N*(Q*)) for n= 1,2. (3.30) 

The compact form of this result makes it straightforward to compute the transport coeffi- 
cients of any lattice gas. Note that knowledge of the conserved quantities and the concomi- 
tant hydrodynamic modes is sufficient to predict the form of this equation. 



3.8 Ordering the Conserved Quantities 



There is one additional technicality that we need to discuss before presenting examples of this 
formalism. In many situations of interest, the conserved quantities themselves are ordered 
in the expansion parameter, e. For example, in an incompressible fluid, the hydrodynamic 
density is assumed to vary by C(e 2 ) from a constant background value, and the hydrodynamic 
velocity is assumed to be 0(e) (low Mach number). 
Consider the general ordering, 



(3.31) 



where it is assumed that the zero-order value, Qq, is always independent of position and 
time. We can then expand the Fermi-Dirac equilibrium as follows: 



«)+E 



dm 



2 ,Uh dQ^dQ" 



(en 
1 



00 v 



K)4g^\ v QtQ\ 



+^ E Ko^-Ko)a-^omh^QiQi 

where we have defined the lowest-order Chapman-Enskog equilibrium, 

and we note that the Fermi metrics and connections in the above expression are now defined 
in terms of Nq , rather than Nq. 
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To make it possible to incorporate this sort of ordering into the above formalism, we 
return to the Chapman- Enskog solution for the Boltzmann distribution, Eq. (|3.25 ), insert 
the above expansion of Nq(Q*), multiply by the fc-fold outer product of the e* vectors and 
contract with gf . In the interests of simplifying the algebra a bit, we shall assume that 
C\ = in this subsection; we do, however, continue to include C2 in our analysis. The 
assumption that G\ — is not in any way essential, but it serves to keep the algebra in 
check. We get 



N l (Q*) 





This result can be directly inserted in the conservation equation, Eq. (|3.2|) 
to 0(e 2 ). After cancelling a factor of e 2 throughout the equation, we get 



where we have defined the advection coefficient, 



V ■ (T>^ e ■ VQf ) + for fieH, 



retaining terms 



(3.32) 



E 



Q1Q1 



(3.33) 



where the diffusion coefficient and source term are still given by Eqs. ( |3.28[) and (|3.29| ), 
respectively, and where all quantities are now understood to be evaluated at A^ . Note that 
the diffusion coefficient is now a function only of Qq, and hence strictly independent of space 
and time. The advection coefficient, on the other hand, is generally quadratic in Qi and 
linear in Q 2 . 

Also note that the advection operator has an 0(1/ e) term and an 0(1) term. If the 
0(1/ e) term does not vanish, it is the dominant term in the equation. In this situation, the 
hydrodynamic equation reduces to the zero-order conservation equation, 



0. 



(3.34) 



Finally, note that knowledge of the conservation laws and concomitant hydrodynamic 
modes is sufficient to predict the form of this equation and compute the advection coefficient. 
Only the computation of the diffusion coefficient requires knowledge of the kinetic modes in 
this ordering scheme. 
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4 Examples of the Boltzmann / 
Chap man- Enskog Analysis 

In this section, we present four examples of the Boltzmann / Chapman-Enskog formalism de- 
scribed in the previous two sections. In each case, we work out the form of the hydrodynamic 
equation, and the transport coefficient (s) predicted by the theory. 

The first example is a one-dimensional, stochastic, diffusive lattice gas (1D3P). The 
second is a lattice gas for Burgers' equation || in one dimension. The third is a two- 
dimensional lattice gas (2D4P) giving rise to coupled diffusion equations. Finally, we consider 
lattice fluid models; we consider a general class of such models, and describe in detail a 
particular model, the FHP-I lattice gas fluid. 

Later in this paper, after we present the complete lattice gas kinetic theory, we shall 
return to these examples, discuss experimentally observed discrepencies from the Boltzmann 
theory, and show how they are explained by the new theory. As we shall see, each example 
has its own unique and interesting features in this regard. 

4.1 The 1D3P Lattice Gas 

As a first example, we consider a diffusive lattice gas model in one dimension (D = 1). The 
model has three bits per site (n = 3), corresponding to the presence or absence of left-moving, 
stationary, and right-moving particles, respectively. These bits are denoted by the respective 
elements of the set B = { — , 0, +}. Collisions occur only if exactly two particles enter a site. 
If we denote the two-particle states by + = {— , 0}, = {— , +}, and + = {0, — }, then the 
nontrivial elements of the state transition table can be written 



a( 


s -> s') 


s' 






+ 









+ 


1 -n p 


n p {l - n r ) 


n p n r 


s 





n p n r 


l-n p 


n p {l - n r ) 






n p (l - n r ) 


n p n r 


l-n p 



The bits n p and n r are random bitsQ which are sampled separately at each lattice site and at 
each timestep with average values (n p ) = 2p and (n r ) = 1/2. Here, the parameter p 6 [0, |] 
may be thought of as the probability of collision from, e.g., + to 0. The value of the bit n p 
effectively determines whether or not a collision will occur, and that of n r determines which 
of the two possible outcomes will result. 

Note that these collisions obey semi-detailed balance, since the columns of the above table 
sum to unity. Also note that they conserve particles (n c = 1); it is the particle density that 
will obey the macroscopic diffusion equation. The coefficients q\ for the conserved quantity 
are 

qi = ql = q\ = 1- (4.1) 

8 Note that r and p are not indices here, but simply labels for the random bits. 



23 



The collision operator is given by Eq. 



cHn*) = n p \n i n i+l n i+2 - n r nV+ 1 n i+2 - r^W+^+si 



n n n n L 



(4.2) 



for each % G B, where n l = 1 — n l denotes the complement of a bit, and where the addition 
of integers to i is understood to increment i through the set B in cyclic fashion. According 
to Eq. we observe that 

= C " + c° + c + . 

We now consider the ensemble average of this system. The ensemble-averaged state 
transition table is 



A(s -> s') 


s' 






+ 









+ 


l-2p 


V 


p 


s 





V 


l-2p 


p 






V 


p 


1 - 2p 



and the ensemble-averaged collision operator in the Boltzmann approximation is given by 
Eq. ( |2.10D which reads 



ji Ari+2 



C l (N*) = 2pN l+l N l+I -pN l N l+L -pN l N 

We do not order the collision operator in this lattice gas, so that C % = Cq. Using Eq. ( |4.1|) , 
we calculate the local Fermi-Dirac equilibrium from Eq. (|3.7f) , 



No = K = iV + 



1 + e- 



(4.3) 



where a parametrizes the distribution function, and / is the mean occupation number. Note 
that the total density is given in terms of a (and /) by 



Q 



l + e- 



3/, 



which is the analog of Eq. ( |3.9|) . 

We now perform the Chapman-Enskog perturbative analysis. The Jacobian of the colli- 
sion operator at equilibrium is 



J 



This matrix has eigenvalues 
and 



' ~2pf Pf Pf 
Pf ~2pf Pf 
V Pf Pf ~2pf 



A 2 = A 3 = -3p/, 
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with corresponding left eigenvectors 



-1 +1 +1 
-1 +1 
-1 +2 -1 



(4.4) 



and right eigenvectors 



Qi 



q-2 







f 


-1 








+2 








-1 



It is easy to check that these eigenvectors satisfy Eq. ( |3.22| ). 

We may now construct the Fermi metric. There is only one hydrodynamic mode, so g^ u 
is a 1-by-l matrix. From Eq. Q3.ll ) we see that its element is 



U = 3/(1-/), 



and hence 



9n 



1 



3/(1-/)' 

To get the generalized Fermi metric, we write c* = ce* where 

e = —x, e° = 0, e + = +x, 

with x being a unit vector, and c the lattice spacing. Using Eq. ( |4.4|) in Eq. (|3.15 ), it follows 
that 



;i) n = o, g (i) 21 = 2/(i-/)x, g (i) 31 = o, 



and 



g(2) n = +2/(1 - /)xx, g(2) 21 = 0, g(2) 31 = -2/(1 - /)xx. 
Using gn to lower indices, we find 

g(l) 1 1 = 0, g(l) 2 1 = fx, g(l)\ = 

and 

g(2) 1 1 = +|xx, g(2) 2 1 = 0, g(2) 3 1 = -|xx. 
Finally, the generalized Kronecker delta, 5(1)^, given by Eq. ( p.23|) , has components 

6{1)\ = ±, *(1) 1 8 = 0. 



The conserved quantity, Q, is not ordered, so we can use the results of Subsection pTF . 
The collision operator was not ordered, so C\ = = 0; Eq. ( |3.26| ) thus tells us that there is 
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no advection or source term in the final hydrodynamic equation. The diffusivity is given by 
Eq. (gjp , 



1 At 



2. ( _ A .) 2 glJ 1 



c 



LuGK 

2 , 2 



3 At V-A 2 
Pxx. 



- 1 XX 



Writing V = xJ-, Eq. ( |3.26|) gives us the hydrodynamic equation, 



dt dx \ dx . 



where the scalar diffusivity is 



V = 



,2 , 2 



3At V-A 2 



: 2 / 2 



3 At V 3pf 



1 . 



Since p lies in [0, ~] and / lies in [0, 1], it follows that V is always positive. 
4.2 The Burgers' Equation Lattice Gas 

The next example is a lattice gas for the one-dimensional (D — 1) Burgers' equation [|J. 
This model has two bits per site (n = 2) corresponding to particles moving left and right. 
The bits are denoted by the elements of the set B = { — ,+}. Collisions occur only when 
exactly one particle enters a site, from either direction. The result of a collision is then one 
particle leaving to the left (state { — }) with probability (1 — a)/2, or to the right (state {+}) 
with probability (1 + a)/2, regardless of the direction of the incoming particle. 
The state transition table is thus given by 



a(s — > s') 


s' 




+ 


s 




1 -n r 


n r 


+ 


l-n r 


n r 



where n r is a random bit with mean (n r ) = (1 + a)/2. Note that particles are conserved, so 

q- = q l + = i. 



The collision operator is given by Eq. 

^{n*) = ± 
We observe that, in accordance with Eq. (P2T 

= cT + c 

The ensemble-averaged state transition table is 



n r n + n — n r n + n 
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A(s -> s') 


s' 




+ 


s 




1-a 
2 


I+a 
2 


+ 


1-a 
2 


1+a 

2 



and the ensemble-averaged collision operator in the Boltzmann approximation is given by 
Eq. ( |2.10D , which reads 



C±(jV* 



± 



- jV+)AT - i— -JV+(1 - N- 



Note that the collisions do not satisfy semi-detailed balance unless a = 0. We therefore 
restrict attention to small values of the bias (a ~ O(e)), and order the Boltzmann collision 
operator as follows: 

C±(N*) = ±~(N~-N + ), 

Ct(N*) = ±^{N- + N + ) =f aN~N + . 

Note that both Cq and conserve particles exactly, and that Cq satisfies semi-detailed 
balance. 

The local Fermi-Dirac equilibrium is given by 



1 + e- 



where a parametrizes the distribution function, and / is the mean occupation number. The 
total density is given in terms of a (and /) by 



Q 



2/- 



1 + e~ a 

The Jacobian of the collision operator at equilibrium is 



J 



This matrix has eigenvalues 



with corresponding left eigenvectors 



2 ~ 2 
1 _ 1 

'2 2 



A 1 =0, A 2 = -1, 

+1 +1 
+1 -1 



and right eigenvectors 



ft = if t\ 



i [ +i 

92= 2 



1 
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Once again, it is easy to check that these eigenvectors satisfy Eq. 

Again, there is only one hydrodynamic mode. The Fermi metric components are given 

by 



9 U 
9n 



i) n =o, g (ir = 2/(1 - /)* 



2/(1-/), 

1 

2/(1-/)' 
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(2) 11 = 2/(1 - /)xx, g(2) 



:i) 1 i=0, g(l)\ 
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0. 



g(2) 1 1 =xx, g(2) 2 1 = 0, 
and the generalized Kronecker delta, 5(1)^, can again be calculated, 

S(l)\ = x. 

Once again, the conserved quantity, Q, is not ordered, so we can use the results of 
Subsection \i.7\ . This time, however, the collision operator is ordered. The first-order collision 
operator respects the conserved quantity, as required, so C\ =0. The component Cf, on the 
other hand is nonzero, so there will be an advective term in the hydrodynamic equation. 
From Eq. ( |3.30| ), we find 



a{N+ + Nq - 2N+N~0) 
2a/(l - /) 

aQ(l- ^ 



The advection coefficient is then given by Eq. ( p. 27] ), 



At^ V 2 



—A 1 



and the diffusivity is given by Eq. (|3.28|) , 



c 

At 

<? 

2At 



v 6(l)\®g(l)\ 1 i 



XX. 
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Writing V = x^, Eq. ( p. 26 ) gives us the hydro dynamic equation, 



3Q + dA = d_ fdQ" 
dt dx dx \ dx , 



where the scalar advection coefficient is 



and the scalar diffusivity is 



v 



2At 



Finally, if we make the change of variables, u = ^(1 — Q), this becomes Burgers' equation 
in the more familiar form 

du ^ du j)d 2u 
dt dx dx 2 



4.3 The 2D4P Lattice Gas 



Next, we consider a diffusive lattice gas model [30] in two spatial dimensions (D = 2). The 
lattice for this model is the standard two-dimensional cartesian lattice. The lattice gas has 
four bits per site (n = 4), corresponding to the presence or absence of particles moving 
along each of the four unit vectors in the lattice, which we will refer to as east, north, west, 
and south. These bits are denoted by the elements of the set B = {1,2,3,4}, respectively. 
Collisions occur if and only if exactly two particles enter a site at right angles, and are 
effected by taking the complement of all four bits at such a site. For convenience, we will use 
a diagrammatic notation for the set of states 5* = {s C B}. In this diagrammatic notation, 
each nonzero bit in a state s is denoted by a line segment emanating from a common vertex. 
Thus, for example, the state s = {2, 3} would be denoted by the symbol 4 . Using this 
notation, the nontrivial elements of the state transition table can be written 



a( 


s -> s') 




s 










i- 


4 


-t 


r 




i- 








1 





s 


4 











1 




"T 


1 













r 





1 









This lattice gas obeys semi-detailed balance, since the columns of the above collision table 
sum to unity. In fact, the lattice gas obeys detailed balance, since the transition matrix is 
symmetric. This collision rule separately conserves the number of east /west-moving particles 
and the number of north/south- moving particles (n c = 2). Thus, the coefficients ql and q 2 
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for the conserved quantities are given by 



Si 1 


= il 


= i 


«l 


= d 


= 


4 


= 4 


= i 


«1 


= <& 


= 



The collision operator is given by Eq. (|2.4| 



(4.5) 



-n i n i+1 n i+2 n i+3 - n l n i+l n i+2 n i+z 
= (y +2 _ n i) (y +1 + n i+3 - n i+1 n i+3 ) , 

where i G -B and the addition of integers to % is taken modulo 4. Since the lattice gas is 
deterministic, we have A(s — > s') = a(s — > s'), and the ensemble-averaged collision operator 
in the Boltzmann approximation is given by C l = d; that is 

C { {N*) = C l {N*) = (iV i+2 - JSr) (iV m + N i+3 - N i+1 N i+3 ) . 
Using Eq. (O), we calculate the local Fermi-Dirac equilibrium from Eq. ( ft.7p , 



o u 1 + e~ ai 

n2 = n2 = — - — = u, 

o o i + e -«a 

where a% and a 2 , or equivalently \x and z/, parametrize the equilibrium distribution function. 
The east /west and north/south densities are given by 

1 2 

Q = z = 2fi 

1 + e" ai 

Q 2 = = 2i/. 

We now perform the Chapman-Enskog perturbative analysis. The Jacobian of the colli- 
sion operator at equilibrium is 



/ -A(v) +A(u) \ 

-AO) +A(/i) 
+A(v) -A(y) 

V +A(ji) -A(/x) / 

where 

AO) = 2z(l - z). 
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This matrix has eigenvalues 



A 1 
A 3 
A 4 



A 2 = 

-2A(/i) 
-2A(z/), 



with corresponding left eigenvectors 



q l = [ 


' +1 


+1 





q 2 =[ 


) +1 





+1 


q 3 = 


} +1 


-1 





q'={ 


> +1 





-1 



(4.6) 



and right eigenvectors 









( °\ 








( °\ 









+ 1 









+ 1 


qi = \ 


+ 1 


qi = \ 





93 = \ 


-1 


q±=\ 







V o J 




V +1 J 




V o J 




V -i / 



We now construct the Fermi metric. There are two hydrodynamic modes, 
represented by a 2-by-2 matrix. From Eq. (|3.11| ) we see that 



so g^ v is 



and hence 



fJ,U 



9/w 



A(/i) 
K{v) 



1 

AM 





To get the generalized Fermi metric we write c* = ce* with 



+x 
—x 



e z = +y 
e 4 = -y, 



where x and y are the unit vectors in the x and y directions, respectively. Using 
in Eq. ( |3.15| ), it follows that the useful nonvanishing components of g(l) are 



Eq. O 



5(1) 



while those of g(2) are 



1(2) 

;(2) 



ii 



22 



A(/i)x 
A(i/)y, 



A(yu)xx 
A(^)yy. 
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Finally, the generalized Kronecker delta, 5(1)^, given by Eq. ( |3.23 ), has components 



The conserved quantities are not ordered, so we can use the results of Subsection 
The collision operator is not ordered for this lattice gas, so C{ = 0{ = and there are 
no advection or source terms in the hydrodynamic equations. The diffusivity is given by 
Eq. O, 



T>» = -— 



c 

At 



6(1) 



( V(v, /i)xx 

I 

where we have defined the scalar diffusivities, 

c 2 / 2 



;-A") 





1 

2 ! 



and 



V(ji, v) 



2At 



2At 



-A3) 



1 



2At \2fi(l-n) 



•A 4 ) 



2At \2u(l 



Writing V = + y^, Eq. (|3.26| ) gives us the pair of hydrodynamic equations, 

dfi d ( dfj,\ 

dv d / . -,dv\ 

Since \x and v lie in [0, 1], it follows that the diffusivities are always positive. 

Note that the diffusivity v) of north/south particles depends only on the density \i of 
east /west particles, and vice versa. In the context of the Boltzmann theory, this is because 
north/south particles scatter only from east/west particles, and vice versa. In Section ||] 
we compute the effects of correlations on these diffusivities, and find that both renormalized 
diffusivities depend upon both particle densities; this is the reason for including the functional 
dependence on both densities in T>. 



4.4 Lattice Gas Fluids 

Finally, we consider a class of lattice gases that have been widely used in recent years for the 
simulation of incompressible Navier-Stokes fluids [1], 0]. Such models exist in dimensions two, 
three, and higher. They conserve mass and momentum, as is appropriate for the Navier- 
Stokes equations. 
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Figure 1: FHP-I Collision Rules 

The nontrivial collision rules for one of the simplest models of this type, known as the 
FHP-I model, are shown in Fig. |l|. The FHP-I model is defined in two dimensions on a 
hexagonal lattice. Note that the two-body collisions have two possible outcomes. We may 
choose between these by a random bit, n r , with mean |. Variants of this model exist with 
other three-body and four-body collisions, with rest particles, etc. We use the above model 
because it is simple but nevertheless fully illustrative for our purposes. Inclusion of the 
three-body symmetric collision is essential because there would otherwise be three conserved 
components of momentum. 

We will now describe the general class of lattice gas fluids and carry out the Chapman- 
Enskog analysis for systems of this type. Following the general analysis we will return to 
the specific example of the FHP-I lattice gas. We consider a general lattice gas with a 
collision rule preserving particle number and momentum in each direction on a lattice in D 
dimensions. We assume that bit i at a site represents the presence or absence of a particle 
of unit massf] and momentum c l /Ai. We can then write the (ensemble-averaged) mass and 
momentum densities as 

n 
i=l 

and 

1=1 AAt 

respectively. The full set of ensemble-averaged conserved quantities for this problem is thus 




Note that this is a (D + l)-component column vector, since there are D components of 
conserved momentum density and one conserved mass density. We shall abuse notation 

9 The restriction to unit mass particles is made only to simplify this presentation, and is not in any way 
essential. Many lattice gas fluid models allow for particles of different masses. 
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by sometimes using p and u as indices; thus, we separate the above equation into the 
components, Q p = p, and Q u = u. Note that Q u refers to D distinct components of Q. 

We can at once identify the left hydrodynamic eigenvectors of the system. They are 
gf = 1 and qf = ^7. To simplify the presentation, we use natural lattice units (c = At = 1) 
throughout this subsection; we can always reintroduce c and At later by dimensional analysis 
considerations. Thus, we write 



1 and qf = e* 



For an incompressible fluid, the conserved densities are ordered in the expansion param- 
eter e as follows fl3lf , 



p = p + e p 2 



u 



eui. 



The second of these equations states that the Mach number is of the same order as the 
smallness parameter used in the Chapman- Enskog analysis; that is, the Mach number scales 
as the Knudsen number. The first equation says that the density fluctuations are smaller 
still - they go as the Mach number squared. 

The ordering of the conserved quantities, Eq. ( |3.31|) , for this system is then, 



Q 




Thus, the zero-order Fermi-Dirac equilibrium is found by considering only q p in Eq. (|3] 

We get 



1 



00 



1 + e~ a 

where / is the mean occupation number of the zero-order equilibrium. Note that / is strictly 
constant - independent of spatial position and time. 

To proceed, it is necessary to impose some requirements on the lattice. Let us assume 
that the lattice is a regular lattice, and that all tensors formed from outer products of the 



lattice vectors are isotropic through the fourth rank |2T| • That is, we demand that the lattice 
vectors be such that 



n 

i=i 

n 

1=1 

n 

i=i 

n 

]T eVeV 

i=l 







±1 

D 



11 



D(D + 2) 



(4.7) 
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where we have defined 

Qijki = 5ij5ki + SikSji + Sudjk- 

We can now construct the Fermi metric. There are D + 1 hydrodynamic modes, so g^ u 
is a (D + l)-by-(.D + 1) matrix. From Eqs. (|3.11|) and ( [4.7|) , we see that it is 



<T = n/(l-/)| 



The inverse metric is then 



n/(l - /) 



1 
±1 



1 
L>1 



Since the conserved quantities are ordered, we are going to need the Fermi connection as 
well. From Eqs. ( |3.14| ) and (|4.7|) , we see that it is 



;i-2/) 



fJtU 



2n/(l - /) 



1 
Dl 



fJiV 



1 

1 



where we have used a row of matrices to represent the variation of the three indices; we 
shall occasionally use this representation for third-rank quantities when it does not cause 
ambiguity. 

Next, we can compute the components of the generalized Fermi metric with hydrody- 
namic indices. From Eqs. ( |3.15|) and (|3.16|) , we have 



c*r =/(i-/) i 



Ei<8> fc+2 e l 



Likewise, the components of the generalized Fermi connection are 

(1-2/) 



2n 2 /(l - /) 

E* ® fe e* DJ2i® 



fc+i 



E*® fc+1 e* £>Ei<S> 



k+2 



e' 
e' 



V ^E 4 ® fc+ 'e J D 2 Ei 

All of the components of these objects which we will need can be evaluated from Eqs. (|4.7|) . 

Eq. (|3.32| ) now tells us the form that the hydrodynamic equations will take. We first 
examine the equation for conservation of mass. Setting the index \i to p, we see that the 
0(1/ e) advection term survives, and all the 0(1) terms vanish; the next surviving terms are 
0(e). Thus, looking to Eq. ( |3.34| ), we get the equation 



V 



o. 
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or, upon simplification, 

V ■ ui = 0, (4.8) 

to 0(e 2 ). We recognize this as expressing the condition that the velocity field must have 
zero divergence in the incompressible limit. 

Turning attention next to the equation for conservation of momentum, we set the index 
H to u. This time we find that the 0(1/ e) part of the advection term vanishes, so the 
hydrodynamic equation is given by the 0(1) terms in Eq. (|3.32|) . For this situation, Eq. ( |3.33| ) 
for the advection coefficient reads 



= g(i) u p g^+[g(i) u ^ uu + g(i) u u r l 

= Pl+p(/)ui Ul , 



rm u . 



: UiUi. 



where we have defined the factor 



and the pressure 



9(f) 



P 



_ (1-2/) 
nf(l - f) 

p2 g{f)u\ 



(4.9) 



D D + 2 1 
written in lattice units. 

Next we turn our attention to the diffusivity tensor, T> u u , given by Eq. (|3.28| ). We have 



*(l) u „ = q?e>d-- 

g(l)"„ = Puug(l) 

g(2) u u = uu g(2) 
so Eq. (|3.28|) becomes 

V Ul h (D + 2)At\ 



n ' 



D JL 



n 



E 

3=1 



e J e J e J e J 



\D(D + 2) 



n 



erefe^ef 



D + 2 ' 



(4.10) 



where we have made it clear that T> u u is a fourth-rank object by explicity writing the four 
spatial indices, i, j, k, and I, and where there is an implicit sum over m and p. Now, 
since the principal fourth-rank tensor constructed from our lattice vectors (see Eq. (|4.7| )) 
is isotropic, and since we have done nothing else to break isotropy, this fourth-rank tensor 
must be isotropic as well. That means that it must be of the form 



(p Ui u y k = v$u$jk + aSijSik + f35 ik 5, 



(4.11) 
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so that the diffusive term on the right-hand side of the hydrodynamic equation becomes 

VJ-D^ V [y kUl = uV 2 u t + V 4 {(a + (3) V ■ u] . 



Note that the second term on the right vanishes due to the incompressibility condition, 
Eq. ( |4.8|) , and v emerges as the shear viscosity in the hydrodynamic equation, Eq. (|3.32| ), 
which can now be written 



dt 



+ g(f)ui ■ Vui = -VP + z/V u. (4.12) 



With the exception of the pathological factor g(f) in the inertial term, which we will address 
below, we recognize this as the incompressible Navier-Stokes equation. 

To get a closed expression for the viscosity, v, we take the following traces of Eq. ( |4.11| ): 

(-D u \) j . = Da + DP + D 2 u 
(p Ui u )[ = Da + D 2 f3 + Du 
= D 2 a + D(3 + Du. 

These three equations may be solved for v to yield 



D(D-l)(D-2) 

Inserting Eq. ( |4.10| ), after a bit of algebra we get the following result for the viscosity of a 
lattice fluid: 

( D v C(e m -e") 2 g- l\ 
(D + 2)At{n(D-l)^ K (-A") 1)' 1 ' ' 



We can compare this result with that of Henon who writes |2T 

C 2 fy 1 

V ~ (D + 2) At [ Xh ~2 



where the quantity Xh is given by 



E Ms - s')r-\l - fr^- 1 [Y(s) - Y(s f )\ : [Y(s) - Y( S ')] 



Xh 2n 

where Y(s) is the traceless part of X(s 



Y(s) = X(s) - 
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and X(s) is in turn given by 



x(s) = sj&i&i 



with trace 



p = Tr[X( S )]=$>' 



Note that our analysis has yielded an alternative expression for Henon's quantity, 



A 



D 



H 



n(D - I) 



E 



Q'n 



So far, we have assumed only that the lattice vectors satisfy Eqs. ( |4.7| ), and that the 
particles each have unit mass and speed. To evaluate Eq. ( |4.13|) for the shear viscosity, 
however, it is necessary to specify a particular collision rule. 

Thus, we now specialize to the FHP-I lattice gas, in two dimensions (D = 2), with six 
bits per site (n — 6), and with lattice vectors given by 



cos 



2ttj\ 
6 ) 



sin 



2vrj 



We then have 



(e^-e fc ) ; 
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The collision rules are illustrated in Fig. 0. We form the Jacobian of the collision operator 
from Eq. ( ^.19| ). Because the collision rules are invariant under cyclic interchange of the 
lattice vectors, the Jacobian is a circulant matrix, so it is particularly easy to write down its 
eigenvalues and eigenvectors. (For more detail on this point, see |[20|| .) We find 



A 1 


= 




A 2 


= 




A 3 


= 




A 4 


= -3/(1- 


/) 3 


A 5 


= -6/ 2 (l- 


-f) 2 


A 6 


= -3/(1- 


/) 3 , 
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with corresponding left eigenvectors 



01 = 


( + 






02 = 


2 






03 = 


\( 






04 = 


i( 


05 = 


( + 






06 = 


2 



and right eigenvectors 



/ + 1 \ 

+i 

+i 
+i 
+i 

V +i J 

( +2 \ 
-1 
-1 
+2 
-1 

"I / 



1 +1 +1 +1 +1 +1 



0+1+10 -1 -1 



-2 +1 -1 -2 -1 +1 



-2 -1 -1 +2 -1 -1 



1 -1 +1 -1 +1 -1 



0+1-10 +1 -1 



/ o\ 

+1 
+1 


-1 

V -i J 



- if 



6 



/ +1 \ 



„6 _ v^S 



V -i / 



/ +2 \ 

-1 

-1 
-2 
-1 

V +i J 

\ 

-1 
-1 


-1 

-1 / 



V 



It is now a simple matter to plug these quantities into Eq. ([4.13 ) to obtain the shear 
viscosity^ 



1 



1 



1=4- 



1 



8At V(-A 4 ) ' (-A 6 ) 7 At\12/(1-/) 3 8/' 

which is the expected result for the viscosity of the FHP-I lattice gas under the Boltzmann 
approximation |2(J Q . 

To conclude this subsection, we return to consider the pathological factor, g(f), that 
appears in front of the inertial term of the Navier-Stokes equation, Eq. ( |4.12|) . Note that 
the convective derivative operator, d t + u ■ V, is Galilean invariant, since it retains its form 
under a Galilean transformation, x — > x' + Vt' and t — > t'. Thus, the presence of the g(f) 
factor in the inertial term is reflective of a breakdown of Galilean invariance. As has been 
pointed out by numerous authors (see, e.g., [0, |2y, 0), this is not surprising since the lattice 
itself constitutes a preferred Galilean frame of reference. 
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We now restore the dimensioned quantities, c and At. 
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In practical simulations of incompressible fluids, this factor is not a problem, since it is 
constant and can be removed by a simple rescaling of either the dependent or independent 
variables. Similarly, the rather unphysical equation of state, Eq. (|4.9|), is not a problem, since 
the equation of state is irrelevant in the incompressible limit. Efforts to extend lattice gas 
methods to treat compressible fluids, however, must deal with these problems. Techniques 
are known for doing this, but they are outside the scope of this paper. 
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5 Exact Analysis 



We will now proceed to develop an exact description of hydrodynamic behaviour in the scal- 
ing limit, dropping the molecular chaos assumption and including effects due to correlations. 
As mentioned in Section [I], the diagrammatic formalism we develop here is similar in many 
ways to the analogous formalism for continuum kinetic theory. Unlike continuous systems, 
however, the discrete nature of lattice gases allows us to explicitly write the complete set of 
terms which contribute correlations over a finite time interval as a sum over a finite number 
of diagrams. The discretization of lattice gases also changes the nature of the vertices in cor- 
relation diagrams. In a lattice gas system, the vertices represent correlation interactions at a 
single lattice site, and can be simply calculated from the time-development equation. There 
are a finite number of distinct vertex types, corresponding to correlated particles arriving at 
and departing from a single lattice point at a single timestep. 



In Subsection |5.1| we generalize our notation slightly to deal with arbitrary sets of par- 
ticles on the lattice. We begin developing the diagrammatic formalism for lattice gases in 
Subsection |5]2| by discussing several alternate descriptions of correlations in ensembles. In 



Subsection [573] we express a renormalized version of the Jacobian matrix J % - in terms of an 



infinite series. The terms in this series are factorized into independent contributions from 



each lattice site in Subsection 5.4, and written in diagrammatic notation in Subsection 5.5 



In Subsection 5.6 we describe how several useful approximations, such as the ring appproxi- 



mation, can be calculated in our formalism as a sum over a restricted class of diagrams. In 



Subsection |5.71 , we prove that a fairly wide class of lattice gases have the property that the 
only effect of correlations is to modify the eigenvalues of the J-matrix in terms of which the 
transport coefficients are described. Finally, in Subsection [5~8| we describe the effects of the 
higher-order collision terms Ci,C% on the exact hydrodynamic equations. 

In Sections || and [7|, we apply the techniques of this section to the lattice gases described 
in Section The reader may find it helpful to refer to Sections |6] and |7| for concrete examples 
of the formalism while reading this section. 

5.1 General Notation 

In this subsection, we develop a slightly more general system of notation suitable for de- 
scribing the exact dynamics of a lattice gas. This notation is similar to the matrix form of 
notation used in Section 2; however, we now wish to consider the space of all configurations 
of the full system, rather than simply the set of states at a single lattice point. 

Recall that we can refer to an arbitrary bit of the system by an index a G B. We now 
introduce an advection operator, A b c , which acts on the space of bits in the entire system. 
This operator is an iV by N permutation matrix. It connects bit b with bit c if and only if 
the particle represented by bit c moves into bit b during the advection substep. That is, 

A b f 1 if & = a(i(c),x(c) + c< c ) 
1 otherwise. 
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In terms of this operator, the exact dynamical equation for the lattice gas, Eq. ( |2.1| ), may 
be written 

n\t + At) = A\ (n c (t) + c i(c) (n a( *' x(c)) (t))) , (5.1) 
and the lattice Boltzmann equation, Eq. ( |2.11| ), may be written, 

N\t + At) = A\ [N c (t) + C l(c) (jV a( *' x(c)) (t))) . (5.2) 

We denote the ensemble mean of an arbitrary product of the n a 's by 




, a C B. 



Henceforth, we use the Greek letters a,f3, . . . to denote subsets of the set B and the letters 
/x, z/, . . . to denote subsets of B. We will sometimes use a roman index to denote an index 
set with a single element, as in N a = . Additionally, for quantities subscripted or 
superscripted by a single set, we will sometimes replace the set by its elements, as in N abc = 
j\[{a,b,c} _ a nna ] p i n t of notation, an index set with a circumflex is assumed to have at 
least two elements; i.e., |S| > 2. 

Finally, we generalize the advection operator A to be a permutation matrix A a p acting 
on the 2 N dimensional space of subsets of B. For a fixed set of bits j3 = {b±, . . . , b q } C B, 
if we take a = {ai, . . . , a q } to be the set of particles which j3 goes to under advection; i.e., 
dj = a(i(bj),x(bj) + c*^)), then 

/ _ f 1, a' = a 
P ~ \ o, a'^oc 

Thus, for example, A a b is 1 when a = a(z(6),x(fo) + c 1 ^), and otherwise, in agreement with 
the previous notation. 



5.2 Representations of Correlations 

An ensemble is generally defined to be a distribution on the space of possible configurations 
of the entire lattice gas system. In this subsection, we discuss several alternative descriptions 
of the probability distribution describing an ensemble. 

Given a set S of boolean variables S = {n l ,n 2 , . . . , n }, a probability distribution on 
S can be described in several equivalent ways. The most familiar description is given by 
assigning a probability to each possible set of values for the n*'s; i.e., given any set a C 
{1, 2, . . . , A^}, we define the probability that the corresponding set of n's are equal to 1 and 
the rest are to be 

P a = Probability that (n 1 , n 2 , . . . , n N ) = (x Q (l), x°(2), • • • , X a (N)), 

where 

*■<'> = {£ \%l 



42 



Since there are 2 N such subsets a, and since 

j2p a = i 

a 

is the only constraint, the space of probability distributions on these N variables is (2 N — 1)- 
dimensional. (In fact, it is a {2 N — l)-dimensional simplex.) 

An equivalent description of a probability distribution on S can be given by defining the 
means N a for each possible product of elements of S. (Note that (n 1 ) 2 = n % , so that the 
mean of any product of elements of S is equal to N a for some a.) 

In terms of the P a, s, the means can be expressed as 

N a = J2 pl3 ' ( 5 - 3 ) 

This relationship can be inverted to get 

p<*= £(-1)101-1°^. (5.4) 

f3Da 

The space of allowed values for the means is also 2 N — 1 dimensional, since N 9 = 1. These 
two descriptions of a distribution are equivalent, in the sense that the information contained 
in either description is exactly sufficient to completely specify the distribution. In fact, 
Eqs. ( |5.3|) and show that the probabilities and the multipoint means are related by a 
linear transformation. 

Probability distributions in which the n 1 are distributed independently have means given 

by 

N a = H N a . 

add 

The space of independent distributions on N variables is clearly N dimensional, and is 
parameterized by N a , a e {1, 2, . . . , iV}. The Fermi-Dirac equilibrium (|3.7|) is an example 
of an independent distribution. 

A third description of a distribution on S can be given in terms of connected correlation 
functions, or CCF's For each a C {1, 2, ... , iV}, there is a CCF, which we denote r a . 
It is easiest to define the CCF's implicitly by expressing the means in terms of the CCF's 
through the equation 

A^ = f*(r*)= r Ci r C2 ...r <: % (5.5) 

where tt(oc) is the set of all partitions of a into disjoint subsets, d, . . . , ( q . Explicitly, 

N a = r a 

etc... 
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We will refer to N a (T a ) as an n-mean (n-CCF), when |a| — n. 

Observing that the defining equation for an n-mean in terms of CCF's contains on the 
right-hand side exactly one n-CCF, and otherwise only m-CCF's with m < n, we see that 
the above set of equations can be inverted by induction on n, to get a functional relationship 
of the form 

T Q = g a (N*), (5.6) 
where g and / are inverses. Explicitly, we have 

r a = N a 

pafe _ N ab _ jyajyb 

Y abc = jy abc _ ]\f a ]Sf bc _ j\r b iV ac - N c N ab + 2N a N b N c . 
etc... 



Again, the description of a probability distribution in terms of CCF's is completely 
equivalent to the discriptions in terms of N a, s and P a, s. There are 2^ — 1 independent r a 's, 
as T is not defined. Note that the relationship between the r a 's and the N a, s (or the P a 's) 
is nonlinear. 

The main reason that CCF's will be a useful description for us is that in an independent 
distribution, all n-CCF's are 0, for n > 1. Thus, the distance of a distribution from one which 
is independent is measured by the quantities T a (recall that an index set with a circumflex 
(~) is constrained to have more than one element). In the subsequent analysis, we will find 
both the N a and T a notations to be useful, and we will use the functions / and g to move 
between the two descriptions. 



5.3 Exact Dynamics 



We will now proceed to rewrite the exact dynamical equation for the lattice gas, Eq. ( |5.1| ), 
in a form similar to that of the lattice Boltzmann equation. For most of this section we will 
assume that the collision operator respects the conservation laws exactly and obeys semi- 
detailed balance, so that C % = C l Q and C\ = C\ = 0; in Subsection |5l| we will discuss the 
effect of correlations when we include nonzero C\ and C\. 

Recall that the collision operator in the lattice Boltzmann equation ( |5.2| ) can be linearized, 
as in Eq. ( |3.18|) , to give 



N b (t + At) = A b c (N c (t) + ef^Nl 



(t) 



(5.7) 



In order to describe the macroscopic behavior of the system we need only include the effect 
of the collision operator Cq up to the order e term associated with the Jacobian, since as we 
have seen in Section ^ the higher order terms associated with the collision operator have no 
effect upon the hydrodynamic equations of the system. In fact, the hydrodynamic equations 
derived in Section [| and the associated advective and diffusive transport coefficients, depend 
upon only the Jacobian matrix J^, through its eigenvalues and eigenvectors. What we now 
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wish to show is that if we drop the Boltzmann molecular chaos assumption and analyze the 
exact ensemble-averaged equation of motion for a lattice gas, in the scaling limit we will get 
an equation identical to Eq. (|5.7|) in form, but with a renormalized J-matrix. The exact 
transport coefficients can then be expressed in terms of the eigenvalues and eigenvectors of 
the renormalized J-matrix using precisely the same expressions as in the Boltzmann analysis. 
Furthermore, we find that for a large class of lattice gases, the eigenvectors of the renormal- 
ized J-matrix are unchanged; only the eigenvalues of the matrix undergo renormalization 
due to correlations. 

We begin with the exact time-development equation, Eq. (|5.1|) . By taking the ensemble 
average of the product of this equation over all a in an arbitrary set a C £>, we can write 
the exact equation for an arbitrary multipoint mean at time t + At in terms of multipoint 
means at time t. We have 

N a (t + At) = (l[n a (t + At)) =J2A a Jl[ [n b (t) + S b \n*(x(b),t))}\. (5.8) 

\a£a I 13 \be/3 I 

To express the right-hand side in terms of multipoint means, it will be convenient to rewrite 
this equation in a more compact notation. For a set (3 C £>, let us define Lp to be the subset 
of points in L which contain at least one particle in the set /3; that is, 



{y G L : x(6) = y for some b G f3}. 



Similarly, we define /3 X to be the set of z's corresponding to the particles in f3 at the point x; 
that is, 

Ac = {i e B : a(i,x) G 0}. 

We can now factorize the product appearing on the right-hand side of Eq. Q5.8Q into contri- 
butions from each of the points in Lp, by writing 

n[n b (t)+c^(n*(x(6),t))] = J] II *) + c V(x, *))" ■ 

be/3 xeL/3 «e/3 x 

The innermost product on the right now depends only on quantities at a single site, x. 
The functions d(n*) can be expressed as polynomials in the n*'s of the form 

C >*) = £a4IK 

uCB jev 

where the k l v are coefficients which may depend only on random bits at each lattice site, 
and which are constant for deterministic lattice gases. Thus, we can write 

nb i + c V)]=<ii" i > 

where the quantities v^ v may contain random bits at each site. Taking the ensemble average 
over any such random bits, we get the mean vertex coefficients V^ u 

n = «>- (5.9) 
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The state transition probabilities A(s — > s') may be interpreted as elements of a collision 
matrix on the space of probabilities, P s , in the sense that the post-collision probability of a 
state s' is given by 

Y,Ms^s')P s - (5-10) 

s 

Similarly, the matrix V^ v can be interpreted as a collision matrix on the space of means. 
Using Eqs. Q5.3Q , Q5.4j) , and ( |5.10| ), the matrix V can be related to A through the equation 

where we have identified the state s with the set of bits which are 1 in that state (s C B). 
Clearly, the 2 2n matrix elements V^ v depend only on the sets fi and z/, and on the form of 
the collision operator. In particular, they do not depend on x, or on the values of the n a 's. 
Note that V\ = 5® u , regardless of the specific lattice gas or collision rule. 
Eq. ( |5.8| ) can now be rewritten in the form 

N a (t + At) = A a p K^N^(t), (5.12) 

where K is an operator expressing the complete collision part of the time development, given 
by 

< = n y\- (5-i3) 

xeL/3 

We shall now transform the exact equation of motion, Eq. (|5.12| ), into an equation of motion 
for the CCF's. Using the functions / and g from Eqs. ( |5.5|) and ( |5.6|) to convert from means 
to CCF's and back, Eq. ( |5.12 ) can be rewritten as 



T a (t + At)=g a (A*pK (3 i r(T*)). 

However, note that from the definitions of / and g, a permutation on the bit labels can be 
performed before or after calculating means from CCF's or vice versa, without changing the 
result. Thus, g and A commute, and this equation can be rewritten as 

r°(t + At) =A a ^ /3 (T*), (5.14) 

where 

$^(r*) = gP(K\r(T*)). 

Since we have chosen to expand around an equilibrium which is an independent distribu- 
tion, all of the quantities T a are of order e (recall \a\ > 2). To perform a complete analysis of 
the system in the scaling limit, we need only keep terms of first order in e in these quantities, 
and hence in the expression $^(r*). This is essentially because the conservation equation 
(p.2|) is linear, and is unchanged by the inclusion of correlations; the effect of correlations 
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appears only at order e in Eq. ( p.lp f"!. We can therefore linearize Eq. ( 5.14 ) in an analogous 
fashion to the linearization of the Boltzmann equation in Eq. (|5.7|), to get 



T«(t + At) = JC~ (etfjsrftt) + K?^ft{t) 



where 



0N a 



K 



T dT"< 



(5.15) 



(5.16) 



Similarly, when we include the effects of correlations to order e in Eq. (|5.14 ) for a = {a}, 
the dynamical equation for iV a = T a becomes 



N a (t + At) =A\{N b (t) + e 



1C\ - 5 b c 



Nt(t)+K*^(t)} 



(5.17) 



Note that if we set T a = in this equation, we get back the linearized Boltmann equation 
since 



K b „ - 5 b 



8 



x(6) 



4(6) 



*{c) u 4(c)' (5.18) 

Inserting Eq. (|5.15 ), we can now write the exact equation of motion for the quantities iV a in 
the form of an infinite series, 



N a (t + At) = A\ (V(t) + eJ b c N?(t) 



where 



x(6) 



x(c) 



J 



4(b) 



4(c) + ^\^ff\ + ^ff\^'^ 



+ ...)). 



(5.19) 



(5.20) 



We now have an expression for the mean occupation number of a certain bit of the system 
at position x and time t + At, written as an infinite sum of terms, each of which is a function 
of the quantities N and at nearby lattice sites x' and at previous time steps t'. As we 
consider terms in this series with more and more factors of A1C, the positions and times at 
which these quantities are evaluated will differ from x and t by greater amounts. However, 



for any given term in the series, the means N^ iy \t') and iVr" '(f) can be replaced by 



r a(i,x'). 



N^(t) and iVr' xj (t), and the expression ( j5.20|) will only change by a quantity of order e, 
since spatial derivatives are ordered as e, and temporal derivatives are ordered as e 2 . Such 
a modification for a finite number of terms does not change the behavior of the system in 
the hydrodynamic limit. In fact, it follows that whenever the sum of terms in Eq. ( p.20| ) 
converges on a scale which goes to zero in the hydrodynamic limit, we can drop all the spatial 
and temporal variations in the single-particle means. Thus, Eq. ( |5.19| ) can be rewritten in 
precisely the form of Eq. (|5.7|) , where J is replaced by the renormalized matrix 



a(i,x) 



(5.21) 



11 The ordering of T a at 0(e) will be assumed throughout this paper. Of course, there will be C(e 2 ) 
contributions as well. These terms are irrelevant in the scaling limit unless there is a divergence in their 
coefficients. We know of no lattice gases which have such a divergence at C(e 2 ) that do not already diverge 
at 0(e). 



47 



with all /C's in J evaluated at the point x and the time t. 

It is important to note that the above argument breaks down when the sum ( 5.20 ) is 
divergent. In this case, the effects of large scale variations in the TV^'^'s must be considered. 
In general, for lattice gases where the sum (|5.20| ) is divergent, one must be quite careful 



about the analysis. For certain lattice gases, however, particularly systems in which the 
conserved quantities are ordered, we are interested in expanding around an equilibrium 
which is spatially invariant (for example, the FHP-I lattice gas). In this case, the zero-order 
means can be replaced by their universal values; however, one must still treat the spatial 
variation of the first-order means carefully. 

Now that we have rewritten the exact dynamical equation in a form commensurate with 
the original form of the lattice Boltzmann equation, the renormalized transport coefficients 
for the theory can be related to the eigenvalues of the matrix J 1 - in the same way that the 
original (Boltzmann) transport coefficients were related to the eigenvalues of the matrix J*.. 
Thus, if we can compute the matrix exactly, we can also compute the exact renormalized 
transport coefficients. Most of the rest of this paper is devoted to methods of calculating 
and approximating the matrix J*., and applications to specific lattice gases. 

5.4 Factorization of Collision Operator 

In the next two subsections we will show that expression ( |5.20| ) for J can be written in a 
diagrammatic notation, allowing us to perform a perturbative calculation of J by summing 
over "Feynman diagram"-like objects, where the contribution from each diagram is just the 
product of factors associated with its vertices. The principal observation which allows this 
reduction is the fact that the collision operator fC a g is factorizable. We devote this subsection 
to demonstrating the exact form of this factorization. 

Theorem 2 For fixed a and f3, K, a p can be broken down into a product of contributions from 
distinct vertices. Explicitly, 

^% = II (5-22) 
where the correlation vertex coefficients (CVC's), V, are defined by 

V% X =E("1) KVI f II ) ( II Ni)v>, (5.23) 

w Ve(axV) / \ie(^\/3x) / 

with the sum taken over fi C a x and v D /3 X , v C B. 

Proof: To prove the theorem, we need only compute the derivatives of the means and 
CCF's with repect to one another, evaluated at the equilibrium point. From the definitions 
and (pD, it is clear that when (3 % a, dg a /dN^ = df a /dV p = 0. When (3 C a, it is 



also fairly straightforward to compute from Eq. ( |5.5| ). 



OTP 



n N o a) = n n k (5.24) 

a£(a\/3) xGL Q ie(a x \/3x) 
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We claim that similarly, when (3 C a, the derivatives dg a /dN 13 are given by 
Lemma 1 

= n (-k w ) = n n hvs)- 

a£a\/3 x£i a ie« x \/3x 



Proof of Lemma: From Eqs. ( |5.6|) and (|5.5| ), one has 



dg c 



on? 



a\7 



7Ca\/3 «G7 



Applying this equation repeatedly, we arrive at an expression of the form 

dg c 



on? 



K\a\(5\)\ J] K 

\iea\/3 



where the function h(n) satisfies h(0) = 1 and 



n 



Kn) = -Y.l =1 ^' m jh(n-m), n^O. 
Assuming that h(n') = (— l) n for n' < n, we have 

n— 1 

h(n) = - £ 



fc=0 



n 
A- 



-1 



-1)", 



so by induction, we have proven the lemma. □ 

Substituting Eqs. QpTTBp , and Q5"25p into Eq. flSTTBD , we get 



En n ( 



n v 



7x 



n n ^ 

V xeL c ie(Cx\/3 x ) 



(5.25) 



Since 7^0, clearly L 7 C L a . The fact that = when v 7^ implies that all terms with 
L( ^ L 7 vanish, so we can restrict the sum over ( to only those ( with L( C L„. In this case 
for x ^ L 7 , clearly V 1 ^ = V \ = 1 or V Jx ^ = 0, and for x ^ clearly ( x \ j3 x = 0, so we 
can replace all the products with products over x 6 L a , giving 



= e n n c-^o)v 

7,C xeLc ie(a x \7x) 



n ^ 



o- 



(5.26) 



ie(Cx\/3x) 



where the sum is over all 7 C a and all ( satisfying C L a and £ D (3. For each x, however, 
this means that 7 X and Cx are summed over all 7 X C a x and C x D /5 X . Since the 7 X and Cx 
are independent for different x, a distributive rule can be applied to Eq. (|5.26|) , giving /C°^ 
in exactly the form stated in Theorem 0, so the proof is complete. □ 



49 



5.5 Diagrammatics 

Using the result from the previous subsection, it is possible to express every term in J in 
diagrammatic form. A generic term in J a h is of the form 

A p al JC a p (5.27) 

with a 1 and /3 l fixed (i.e., not summed over), and \a l \, \(3 l \ > 2, except for the endpoints 
where a k = {a} and [3 1 = {b}. 

We define a diagram T by an integer k(T) , which we refer to as the length of the diagram T, 
and a function a^W) where for each r 6 {0,...,fc},aT(r) C B. Geometrically, we associate 
each a G «t( t ) with a virtual particle (VP) moving from (x(a), r) to (x(a) + c l w, r + 1) on 
the lattice A fc+1 = L x {0, . . . , k + 1}. We refer to ^(r) as the set of outgoing VP's for the 
diagram T. 

It is natural to define a corresponding set of incoming VP's for r > by Pt{t) = {b : 
a(i(6),x(6) — c*( 6 )) G «t(t — 1)}. We also define <7t(t) = |«t(t)| to be the total number 
of outgoing VP's for each value of r. Finally, given a diagram T, we can define a weight 
function 

w(T) = u n v- ( ^ T(r)x , 

xeL l<r<fc(T) 

by taking the product of V over all vertices. 

The term ( |5.27|) can now be represented by the diagram T with q;t(t) = « T , where for 
consistency a is defined to be the unique set with A = 1. When r ^ 0, A T = 1, so 
!3t{t) = f3 T for all r. It follows that the contribution from the term (|5.27|) is exactly given 
by W(T). Thus, we can rewrite expression ( |5.20p for J\ as a sum over diagrams 

oo 

J\ = H E ^( T )> (5-28) 

k=l TeT a b (k) 

where in general we define the set of diagrams T a Jk) by 

T a p {k) = {T : k = k(T), a T (l) > 1, for 1 < / < k, a T (k) = a,(3 T (l) = (3}. 

Note that V® v = 5® v) so any diagram with incoming VP's at (x, t) but no outgoing 
VP's has weight zero. For many lattice gases certain other vertex factors V^ u vanish also; 
diagrams with such vertices can be dropped from the sum (|5.28 ). From Eq. (|5.21|) , J can 
now be written as a sum over diagrams in the same fashion as J , 

oo 

J l M) = H E W(T), (5.29) 

fc=lT6T*(x,ft) 

where the set of diagrams to be summed over is given by 

r;.(x,fc)= U T a ^\(k). 

b:i(b)=j 
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r = 1 




t = 2 



r = 3 



Figure 2: k = 3 Diagrams for a ID Lattice Gas with Three Bits/Site 



We will find it useful later to generalize this set of diagrams to the sets of diagrams 

Tt(x,fc) = {T:k = k{T), a T {l)>l, fori <l< k, |L ar(0) | = 1, L QtW = {x}, 
li(a T (k)) = /I, /i(a T (0)) = v}, 

where we have used the notation 



Example: As a simple example of the diagrammatic notation, consider the allowed di- 
agrams for the 1D3P lattice gas considered in Section 2.5. The complete set of diagrams 
needed to compute the k = 3 correction to J° is the set of diagrams Ti - T 4 represented in 
Fig. along with the diagrams achieved by reflecting Ti and T 2 across the x-axis0- The 
weight of diagram T 2 , for example, is W(T 2 ) = V°~V + _V° V + . We will compute the vertex 
factors V M „ for this lattice gas in Section |6A| , and we will see that the contribution from 
diagram T 3 in fact vanishes. 

5.6 Approximations 

We have so far managed to write the exact formula for the hydrodynamic equations in the 
scaling limit only in terms of an infinite formal series. The natural next question to confront 
is whether this series can be summed. We would like to know whether the series is finite, 
and if we cannot sum the full series, at least we would like to find a set of reasonable 
approximations which we can make to truncate the series to one which is summable. The 
questions of convergence are rather difficult, and we will not address them here in full 
generality; in general, the convergence properties of the series depend on the form of the 
conserved quantities in the system. A variety of methods for performing partial sums of 
infinite series of diagrams while retaining physically important terms have been applied to 
related problems in quantum field theory and quantum many-body theory |33|, |33| . We will 
describe here several particular approximation methods which are useful for the kind of series 
which arise for known lattice gases. 

12 Throughout, we shall depict diagrams for one-dimensional lattice gases with a vertical time axis. 



//({ai, . . . , a,-}) = {i(ai), . . . , 
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5.6.1 Short-r and Small-^ Truncations 



The simplest useful approximations involve truncating the sum ( 5.28Q to a finite number 
of terms by putting an upper bound on either the number of timesteps or the number of 
distinct nontrivial vertices allowed in each diagram. In the first case, the expression for the 
renormalized J-matrix is 

k=l TeT\ (x,fc) 

where the diagrams summed over are the same as those summed in Eq. Q5.29 ). Since for 
each fixed value of k there are a finite number of allowed diagrams, this sum is finite. We 
refer to this approximation as the short-T approximation. In the second case, we allow k 
to be arbitrary, but allow only diagrams where the total number of vertices (x, k') with 
nonempty outgoing sets «t(^') x is less than or equal to some fixed number £. We denote 
the sum restricted to these diagrams by J w r A gain, there are only a finite number of such 
diagrams in this sum, which means that this sum must also be finite. This approximation 
is analogous to the weak-coupling expansions in quantum field theory, although in this case 
the coupling constants V^ v are usually not particularly small. The short-r and small-^ 
truncations give good consecutive approximations for many lattice gases. In either of these 
two approximations, the Boltzmann approximation can be recovered, by taking r = 1 or 
£=1. 



5.6.2 BBGKY Truncations 

Another good class of approximations, in which a reduced but still infinite set of diagrams 
is summed, corresponds to truncations of the BBGKY hierarchy of equations. Such an 
approximation involves neglecting g-CCF's with q > n for some fixed value of n. In our 
diagrammatic formalism, this amounts to restricting the sum to diagrams with crik') < n 
for 1 < k' < k. For example, with n — 2, diagram T 3 of Fig. [| would be neglected. 
Whereas the computational complexity of including all diagrams in the complete sum grows 
exponentially in k, that of the truncated BBGKY approximations grow polynomially, and 
are therefore computationally more tractable. 



5.6.3 The Ring Approximation 

The n = 2 version of the BBGKY approximation is closely related to the ring approximation. 
The ring approximation is made by neglecting interactions between two propagating corre- 
lated quantities except at the initial and final vertices of a diagram. It is generally possible 
to calculate a closed-form expression for the infinite sum of diagrams corresponding to this 
approximation. Furthermore, it is usually fairly easy to calculate the asymptotic form of this 
approximation as k —>■ oo. This calculation often captures the most significant part of the 
long-term renormalization effects. In particular, for certain lattice gases which model two 
dimensional fluid systems, the ring approximation diverges logarithmically in \L\, which is in 
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agreement with predictions from other theoretical frameworks E|, and also with observed 
behaviour [ j35fl . 

In Section 7, we will apply the different approximation methods described here to the 
1D3P lattice gas and compare the results from these approximations with experimental 
results. 

5.7 Eigenvalue Renormalization 

Within the framework of the formalism developed in the previous subsections, we can now 
demonstrate that for a large class of lattice gas models, the effects of correlations are to 
renormalize only the eigenvalues of the J-matrix, and not to change the eigenvectors. This 
result follows from a pair of fairly simple theorems. 

Theorem 3 The finite matrix C • = J j — J l j of corrections to the J-matrix can be 
restricted to be a matrix in the space of kinetic eigenvectors of J; i.e., if v 6 H then 

qUJ^j - J)) = (J (T) ) - J)H = o. 







Proof: From Eqs. ( |5.20| ) and ( |5.22| ), it will clearly suffice to show that q^V* 
for every |iC5 with \p,\ > 2. 

We first demonstrate this sufficient condition for the right hydrodynamic eigenvectors q l y . 
We showed in Section [375] that q l u was a right eigenvector of the J matrix with eigenvalue by 
differentiating the identity C 1 (Nq) = 0, which holds for any equilibrium, with respect to the 
parameters a u of the equilibrium. We can similarly show that q l v is a null right eigenvector of 

V M , by using the stability of the Boltzmann equilibria. The stability of the local Boltzmann 
equilibrium tells us that at each lattice site, 



pCB i£p 

Differentiating this equation with respect to the parameter a v of the local Boltzmann equi- 
librium gives 

o = e ( n -*t J K E( n nSH = 

as desired. 

Next, consider the left hydrodynamic eigenvector q\. From Eq. ( |5.23| ), we have 
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Since q\N l is a conserved quantity for any values of N l , we have 
for arbitrary Nq, so gfVJ- = <zfV~ = for all p, C 5 with |/t| > 2.D 

Theorem 4 // t/iere exists a symmetry E of a lattice gas which can be expressed as a com- 
bination of a permutation E x y on the lattice L fixing a point x and an independent per- 
mutation 5j ■ on i/ie 6it sei _B ; where the zero-order Boltzmann equilibrium at x satisfies 
E l -iVo(x) = A^(x) ; i/ien E commutes with J and ; i.e., 

- fjM&k = sy (T)3 fc (x) - j (t) ;-(x)s\ = o 



Proof: The condition that E is a symmetry of the lattice gas with the given product 
structure asserts that E a b = E ^E* x n,) * s a permutation matrix on all the bits of the 
system which can be extended to a permutation matrix S"^ on subsets of B, satisfying 

(E- 1 ) V^ 7 * = A% 

and 

Since by their definitions, the functions / and g connecting means and CCF's are invariant 
under permutations, it follows from Eq. ( |5.16 ) and the invariance of the local Boltzmann 
equilibrium at x that 

(s-r^ 7 E^ = ic% 

where we assume that all zero-order single particle means are evaluated at x, as in the 
diagrammatic expansion. It then follows from Eqs. ( f5.18D , ( 5.20|) , and the fact that 5 a p 
commutes with E that 

(E- 1 ) V 6 c E c d = J\. 

and 

(e-^me* = j*(x), 

proving the first part of the theorem. From the constraint that E leave the point x fixed, we 
have 



^(x) = E^ a(l ' x) 



a(i,w) 
w 

E(s- 1 ) x y (s- 1 )y a( n ( ,, z) s fc i s z v 

w 

e (s^y^'tcM^ 

z 

(E-^.j^x^, 



so the second assertion of the theorem is proven. □ 
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Corollary 1 If the sequence {J j} is convergent as r — > oo, i/ien Theorems [3| and [| /io/d 
m this limit. That is, 



tii^j - Ji j) = ( J ) ~ J )H = °> for v £ H, 



and 



^jJ k k ~ ^ l jJ k k — 0' 



for E* • as m Theorem [|. 



For a wide class of lattice gases, including all the lattice gases described in Section pTT| , 
these two theorems suffice to demonstrate that the eigenvectors of the J matrix coincide with 
those of J, so that the only effect of correlations is to renormalize the eigenvalues of J. The 
essential point is that we can classify the eigenvectors of J by their transformation properties 
under the group of symmetries of the lattice gas. Theorem § asserts that J can only mix 
eigenvectors with identical symmetry properties. Thus, if no two kinetic eigenvectors of J 
share the same symmetry properties, then J must be diagonal with respect to the basis of 
eigenvectors of J. More generally, if the representation of the group of symmetries on the 
space of kinetic eigenvectors of J breaks up into irreducible representations in such a way 
that no irreducible representation appears more than once, then J is diagonal with respect 
to the eigenvectors of J. 

As a simple example, consider the 1D3P diffusive lattice gas. The right kinetic eigenvec- 
tors of ,/ arc 

</2 = - ! o 

+i 




and 

93 = 

Under the symmetry transformation x — > — x, + <-> — , these eigenvectors transform with 
eigenvalues S 2 2 = — f, S 3 3 = +1. Thus, these eigenvectors cannot be mixed by J, and so 
their eigenvalues are separately renormalized. Explicitly, in matrix notation, with respect to 
the basis we have 



J 



By Theorem [|, J commutes with E, so j 2 3 = j 3 2 = 0. 

For most standard lattice g similar analysis of the symmetry properties of the 

eigenvectors of J shows that no irreducible representation of the symmetry group appears 
more than once, so that the eigenvalues of J are renormalized in a straightforward fashion. 
For those lattice gases where this cannot be shown, it is necessary to repeat the entire 
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Chapman-Enskog analysis using the renormalized J matrix. Note that Theorem ^ implies 
that any eigenvectors of J which lie in the same irreducible representation of the symmetry 
group must have identical eigenvalues in the matrices J and also in J. 

5.8 Renormalization Effects and Higher-Order Collision Opera- 
tors 

We conclude this section with a discussion of the effects of the higher-order collision operators 
eCi and e 2 C2 in the full kinetic theory. We show that the second-order term only appears 
in the source term for the hydrodynamic equation as in Eq. ( 3.29|) , and does not generate 



extra renormalization effects. On the other hand, we find that the first-order collision op- 
erator not only appears in the advection coefficient as in Eq. (|3.27| ), but generates a set of 
additional correlations which modify the advection coefficient by effectively renormalizing 
the components of the first-order collision operator itself. 

We begin by discussing the second-order term e 2 C*2. This part of the collision operator 
only appears in the first-order conservation equation ( |3.26|) . The corrections to this term 



due to correlations are of one higher order in e, and can clearly be neglected in the entire 
analysis. Thus, inclusion of this term only generates the source term in Eq. (|3.26| ) in the 



manner described in Subsection 3.7 



Now, we consider the effects of including a first-order term eC\ in the collision operator. 
Recall that this term is restricted to obey the conservation laws, but is not required to satisfy 
semi-detailed balance. Because this part of the collision operator obeys the conservation laws, 
it does not appear directly in the first-order conservation equation. It does appear in the 
linearized Boltzmann equation (|3.18|) . However, corrections to the term in this equation due 



to correlations are again of higher order in e, so that no change is necessary to this equation 
due to renormalization effects. 

At this point, one might imagine that inclusion of the first-order term in the collision 
operator does not necessitate any further modification to the exact hydrodynamic equation 
in the scaling limit other than the effects described in Section 3. However, this is not the 
case. In fact, the inclusion of this term in the collision operator has a nontrivial effect on 
Eq. ( |5.15D , which describes the propagation of correlated quantities in the system. 

Eq. ( [5.15| ) gives an expression to order e for T a at time t + At in terms of quantities 
at time t. In the derivation of Eq. (|5.15| ), we used the fact that the zero-order means Nq 
describe a local Boltzmann equilibrium which does not generate correlations through the 
collision operator Co- When we include the first-order collision operator Ci, we must include 
the fact that €C[(Nq) need not vanish. Thus, the correct form of Eq. ( pM5| ) in this case is 

V s (t + At) = (eK?Nt{t) + K?J^{t) + elA , (5.30) 
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where we have defined I@ = whenever \L~\ > 1, and 



i?=jm S £ n 




e=0 



/ jet 



(5.31) 



whenever L-g = {y} and (3 y = = \i for some y G L; all means iV in this equation are 
evaluated at the point y. In Eq. ( |5.31| ), we have used the collision operator on means V" 
defined by using the collision operator C % = Cl + tC{. However, the CVC's used in Eq. ( |5.30|) 
should still be evaluated with respect to the zero-order collision operator Cq. 

Combining Eqs. ( |5.30| ) and (|5.17|) , and including the first order collision operator as in 
Eq. ( ETT5D , we get 

N a (t + At) = A\ (N b (t) + eJ b c N{{t) + el b ) , 

where 

X b = C f\N*^c{b))) + K, b z{AyP + A^JC^{Ay ? + . . .)). 

By applying the same analysis used in Section 3, we find that the renormalized hydrodynamic 
equations for the theory are of the same form as Eq. fl3.26| ); however, the inclusion of the 
first-order collision term changes the result for the renormalized advection coefficient to be 



A^Q* 



At 



E 



where the renormalized generalized Kronecker delta function and eigenvalues are defined 
with respect to the renormalized J matrix, and where the renormalized collision operator C{ 
is given by 

Cj'(x) = q v i X< i ^\ 

Just as for the renormalized J matrix, we evaluate all means N appearing in CVC's and 1^ 
at the point x. This simplification depends again upon the convergence of the infinite series 
of terms in C". 

We can express the renormalized collision operator in terms of an infinite diagrammatic 
sum, analogous to the sum (|5.29 ) for the renormalized matrix J. Specifically, we have 



X a = C i{a) + J-, 



W{T)I". 



k = l T&T i[a) (x(a):fc) 



57 



6 Examples of Vertices and Renormalization 



We will now apply the methods of the previous section to compute the vertex factors for the 
example lattice gases described in section 4. We also derive expressions for the renormalized 
transport coefficients for all these lattice gases. 



6.1 1D3P Lattice Gas 
6.1.1 Vertices 

Beginning with either the collision operator (|4.2|) or the state transition table (see Subsec- 
tion |4.1| ) for the 1D3P lattice gas, we can calculate the mean vertex coefficients V^ v using 



Eqs. (|5.9| ) and (|5 . 1 1|) . The nonzero mean vertex coefficients are given by 



' B — V — 1 

V) = 5\ 

Vi = p(3^.-l) 

vf = p + ^.(l-3p). 

Note that since the ensemble-averaged collision operator is invariant under permutations 
(relabeling) on the bits, the mean vertex coefficients also have this symmetry. 

Using the equilibrium value / for the mean occupation numbers Nq (see Eq. ( |4.3j )), the 
expression for the CVC's, Eq. ( |5.23|) , reads 

\> a f3 = (_i)M-M f\ a \+W\-M-\P\y» _ 

For instance, we have 

v+ = f(v\ + v^-v -v- )-f 2 (v% + v~ z + v°~ + v-~) 

= -pf(l-f) 

The remaining nonzero CVC's are given by the equations 

V B B = V\ = l 

v] = <*yw(i -f)-pfO—f) 



VV = Kl-/) + ^(l-3p(l-/)). (6.1) 



Note that the CVC's are also symmetric under an arbitrary permutation on the particle 
labels. The nonvanishing correlation vertex factors are depicted graphically in Fig. § only 
a single vertex is shown in each equivalence class under the permutation symmetry. 
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V° = 1 - 2pf 



} 

V~ =pf 




Figure 3: Vertex Factors for 1D3P Lattice Gas 
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6.1.2 Renormalization of Diffusivity 



As was shown in Section |5.7| , the symmetry of the 1D3P lattice gas under spatial inversion 
is sufficient to ensure that the effect of correlations is simply to renormalize the eigenvalues 
A 2 and A 3 . From the Chapman- Enskog analysis, we know that the diffusivity is given in the 
hydrodynamic limit by 

D = f-t- - l) , (6.2) 
3At V-A 2 / 

where A 2 is the eigenvalue of the vector (—1,0,-1-1) in the matrix J. In this subsection we 
will find the set of diagrams which contribute to the renormalization of A 2 . 

Factoring out the initial and closure vertices, we can write the renormalized matrix J as 

^.(x) = 4 + vvv i ;.y|(x), (6.3) 

where 

oo 

fc=lTeT" M (x,/c) 

From the fact that V B = V B { = for all % 6 B, we see that the only nonzero contributions 
to Eq. (|6.3| ) can come from terms of the form VlY\y k .. From the inversion symmetry, we 

I k j 

observe that _ _ 

y r L. = yt 

-k k 

Thus, Eq. ( |6.3|) can be explicitly rewritten as 



\ +1 + +3 +3 Oj + I 0] 0' 



where 



X%. = V\V.. 

kj k j 



To evaluate A 2 , we need now only calculate the eigenvalues of qf and q{ with respect to the 
matrices X 1 -., for each pair of values for I and k. Evaluating 



kj 



q-V% = -g?V~ = 3p, 

= -v-^ = ^p/(i-/), 
^ = -v°^ = o, 



we have 



A 2 + 2qfX i+ ^qiY\ + 2q*X i ^.<&Y^ (6.4) 

,2 . 



-3pf + 9p 2 f(l-f)(Y\-Y^ 
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Diagrammatically, this equation can be expressed as 



-3pf + 9p 2 /(l - /) 



A 
V 



A 
A 



\ 



(6.5) 



/ 



where the notation in brackets indicates summation of the products of all internal vertex fac- 
tors over all diagrams with the depicted initial and final configurations. Together, Eqs. (|6.2j ) 
and ( |6.5[ ), with vertices given in Fig. ^ constitute an exact expression for the diffusivity of 
the 1D3P lattice gas. 



6.2 Burgers' Equation Lattice Gas 
6.2.1 Vertices 

From either the collision operator or the state transition table for this lattice gas (see Sub- 
section |4.2| ) we can calculate the nonzero mean vertex coefficients, which are given by 

1 



n 


- V - 




= —a 


v- B 


= a 


v\ 


1 + a 
2 




1 — a 
2 



Using the equilibrium value / for the mean occupation numbers Nq , we can calculate the 
correlation vertex coefficients. Recall that we calculate the CVC's using only the zero-order 
mean occupation numbers. The CVC's are given by 



v B 

v*. 

v 3 
v B 
v B 



V 

1 

2 



B 



B 



1 



= 0. 



In particular, note that all the CVC's which modify the number of correlated quantities 
are zero. Thus, in this lattice gas, no 0(e) correlations are generated by gradients in Aq, 
and correlations cannot affect the hydrodynamic equation by influencing the single-particle 
means. 
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Because of the first order collision operator C\ which does not satisfy semi- detailed bal- 
ance, correlations might also be generated by the quantities I M . From Eq. ( [5.31| ), we have 



r-B 



0. 



Thus, we find that for this lattice gas, no correlations are generated to 0(e). Furthermore, 
even if correlations existed, they would not couple back to the hydrodynamic equations, since 
V B = 0. It follows that the standard Chapman-Enskog analysis gives the correct results for 
the transport coefficients. In fact, this result was proven using other methods in [26 . 



6.3 2D4P Lattice Gas 
6.3.1 Vertices 



Using the general formulae, Eqs. ( |5.11 ) and ( |5.23|) , to calculate the correlation vertex coeffi- 
cients V a p for the 2D4P lattice gas, we arrive at the values for the CVC's which are depicted 
in Fig. [|. In this figure we have only included a single example of each equivalence class of 
vertices under the symmetry group generated by reflections across either axis. Furthermore, 
we have only included a single example of each pair of vertices which are related by a 90° 
rotation. The CVC's for such vertices are related by exchanging the two particle densities 
\l <-> v. For example, we have 



l-2v + Iv 1 
1 - 2/i + 2^ 2 . 



(6.6) 
(6.7) 



6.3.2 Renormalization of Eigenvalues 



Applying the theorems of Subsection |5J to the 2D4P lattice gas, we find that the renormal- 
ized J-matrix must be of the form 



J 



V 



A 4 /2 





-A 4 /2 


\ 





A 3 /2 





-A 3 /2 


-A 4 /2 





A 4 /2 








-A 3 /2 





A 3 /2 ) 



The renormalized hydrodynamic equations are 

d 



with 



~dt 
dv 
~dt 



f>( ) ^ 

dx \ ' dx / 

dy I ' dy 



2At 



-A 3 ) 
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1 - 2u(l - v) 



2i/(l - u) 



- u)(l - 2v) 



-u)(l- 2u) 2/i 2 (l — /i) 2 



-2 M 2 (1-^) 2 



1-2;/ 



1 - 2j/ 2/i(l - ;u) - 2i/(l - v) 



2z/(l - i/) + 2/i(l - /i) 1 - 2i/(l - i/) - 2/i(l - n) 



- m)(1 - 2/i) m(1 - - 2//) 



-2 



+2 



1 - 2v 



-(1-2^) 



2/i(l - /j) 1 - 2/i(l - A*) 1 

Figure 4: Vertex Factors for 2D4P Lattice Gas 
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and 



2At 



-A 4 ) 



In order to calculate the renormalized eigenvalue A 4 it is only necessary to calculate the 
component J\ of the renormalized J-matrix. Because the lattice gas is invariant under the 
combination of a 90° rotation and the exchange \i <-> u, the eigenvalue A 3 can be calculated 
from A 4 by the exchange of particle densities. 

Using a similar notation to that developed for the 1D3P lattice gas in Section |6|.1.2, we 
can write the renormalized matrix J in the form of Eq. ( |6.3|) . Using the symmetry properties 
of the lattice gas to prove the equality of quantities Y a ^ which are related by reflections across 
each axis, we can simplify the expression for the shift in A 4 to 



5\ 4 = Su{u- l)(2v- l) 2 
+16^0 - l)(2i/- 1 
Similarly, we have 

5A 3 = 8/i(/i- 1)(2//- 1) 
+16/x(//-l)(2//- 



-Y 



Y\ + Y\ 



Y\ 



l<ov 2 {v- if 



y f -y f 



\<ov 2 (v -\ ?{2v - 1) 



+ 16//(/i - if 



Y 



Y 



Y +Y 



+ 16fi 2 (fx-lf(2fx-l) 



Y 



Y 



These equations describe completely the renormalization of the hydrodynamic equations 
due to correlations. Each term Y a ^ corresponds to a set of diagrams with a specific set of 
outgoing and incoming virtual particles at the initial and final vertices of the diagram. 

A rather dramatic simplification of the eigenvalue renormalization equations occurs when 
the equilibrium particle densities /x and v are equal to 1/2. When v = 1/2, we have 



and when /i = 1/2 similarly 



5X 4 



5A 3 



Y 



f 



Y 



4 



(6i 



(6.9) 



It follows that for these particular values of \x and u, the set of diagrams which give a nonzero 
contribution to the renormalization of the eigenvalues is reduced to only those diagrams which 
have 3 outgoing virtual particles at the initial vertex, and 3 incoming virtual particles at 
the final vertex. Thus, for example, the ring and 2-particle BBGKY approximations for this 
lattice gas vanish at the equilibrium described by fj, = v = 1/2. 

In this particular case, the r = 3 short-r approximation is given by Eqs. ( |6.8|) and ( |6.9|) , 
with the first Y term in each expression vanishing (because it is impossible to connect the 
outgoing and incoming particles in two time steps), and with the second arising from a single 
diagram of weight 1/8. It follows that the corrected eigenvalues are equal to —9/8 in this 
simple approximation. This leads to T>/T> = 7/9 = 0.777. . ., which may be compared with 
the experimental value of "Cexpt/^ 5 ~ 0.71 HI- 



64 



J 



■/(I-/) 2 3/(1 - /) 2 (1 + /) 3/(l-/) 2 (l-3/) -/(l - /) 2 (1 - 2/) 



Figure 5: Propagator Vertex Factors for FHP-I Lattice Gas 



In a future paper [Q, we will discuss in more detail the results of summing various 
subsets of diagrams for this lattice gas and compare the results to experimental data. The 
ring kinetic theory for this lattice gas has been worked out in []37| ; in this reference, a ring-like 



approximation is also used to treat noninteracting three-particle correlations, giving nonzero 
correction in the \i — v — 1/2 case. 

6.4 FHP-I Lattice Gas 
6.4.1 Vertices 



Using the general formulae (|5.11| ) and (|5.23| ) to calculate the correlation vertex coefficients 
V a p for the FHP-I lattice gas results in nearly 300 nonvanishing CVC's that are independent 
in the sense that they are not related by symmetries. While it is a straightforward task for a 
symbolic algebra computer program to compute and work with these quantities, it would not 
be useful to present all the results in this paper. Instead, we present only the initial (one-to- 
many) vertices, the propagator (one-to-one) vertices, and the closure (many-to-one) vertices. 
These vertices are sufficient to compute all the diagrams in the kinetic ring approximation, 
and in the r = 3 short-r approximation. 

6.4.2 Renormalization of Eigenvalues 

The only transport coefficient of the FHP-I lattice gas which undergoes renormalization is the 
shear viscosity v. From the general arguments in Subsection |5.7|, we see that the renormalized 



eigenvalues A 4 and A 6 must be equal, as they lie in a single irreducible representation of the 
symmetry group generated by a 60° rotation. In terms of these renormalized eigenvalues of 
the J matrix, the renormalized shear viscosity v is given by 



c 2 



8AtV(-A 4 ) J 

The renormalized eigenvalue A 4 can be extracted from the J matrix by taking the components 

A 4 = 6(J\ + J\). 
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_-/. _ -V 

^/(l-/) 2 -|/(l"/)(l-3/) -(l-/) 2 (l-2/) -/(I-/) 2 




|/(l-/)(3-B/) i(l-/) 2 (l-2/) -/ 2 (1-/) i/(l-/) 2 -/(l-/)(2-3/) 




/(I-/) £(l-/)(2-5/) -1(1 _/)(!_ 2/) -/(I-/) i(l-/)(2-/) -1 + 3/(1-/) 




-/(I-/) -|(l-/)(l-5/) -1(1-/) 2 - 1(1 -/)(! + 2/) 1-3/(1-/) /(I-/) 




-2 + 3/ -/ j(l-/) -1(1-3/) 1-2/ 




Figure 6: Initial Vertex Factors for FHP-I Lattice Gas 
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Figure 7: Closure Vertex Factors for FHP-I Lattice Gas 
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As in the previous examples, we can express the eigenvalue shift in the form 



6\ 4 = 6 f X 1 ^ + X l \Y"~ 

By collecting coefficients of terms Y ~ which are related by symmetry, the expression for the 
eigenvalue shift reduces to a sum over more than a hundred terms. 

We will not explore the diagrammatic expansion of the FHP-I lattice gas further in this 
paper. However, we have used the general expression for the renormalized shear viscosity 
to calculate several simple renormalization effects. In particular, one finds that in the ring 
approximation, the sum over diagrams diverges logarithmically; this result is well-known 
for incompressible fluids in 2-dimensions, and has been verified using other methods for the 
FHP-I lattice gas @. 



68 



7 Approximations and Numerical Results 
for the 1D3P Lattice Gas 



In this section we discuss in detail the sets of diagrams for the 1D3P lattice gas which 
correspond to the various approximation methods described in Section 0. 7, and perform the 
associated partial diagrammatic sums. The goal of this detailed analysis is to use the simple 
1D3P lattice gas as a test case to study the relative accuracy of the different approximations 
and the relative difficulty of computing these approximations. The results for this lattice 
gas hopefully give a good indication of what approaches will lead to useful results for more 
complicated lattice gases. 

Correlations cause corrections to the Boltzmann diffusivity of the 1D3P lattice gas which 
are as large as 5% for certain values of the particle density / and bounce probability p. 
We describe in this section the results of numerical calculations of the partial diagrammatic 
sums corresponding to the various approximation methods, and compare to empirical results 
from computer simulations of this lattice gas. In Subsection we calculate corrections to 
the diffusivity in the short-r approximation. By graphing these corrections and compar- 
ing to experiment, we see that as r increases these approximations give corrections to the 
Boltzmann approximation which seem to converge to the correct experimental values. The 
convergence is slow, however, and since the calculation of these corrections is computation- 
ally quite expensive, it is difficult to estimate the asymptotic value of the diffusivity to a 
high degree of accuracy using these approximations. In Subsection |7.2| we apply the par- 
tial BBGKY summation prescription and again graph the results compared to experimental 
values. For small values of k the fc-particle BBGKY diagrammatic summation converges 
rapidly, and in most ranges for which we have calculated the results, these approximations 
approach monotonically the experimental values as k increases. In this subsection we also 
include a proof that the 2-particle BBGKY approximation converges for arbitrary values 
of p, f > 0. In Subsection [7]3] we observe that in the vicinity of the density value / = 1, 
the corrections to the Boltzmann approximation can be expanded in a power series in the 
variable e = 1 — /. The diagrams which contribute corrections of order e n are the diagrams 
which have n or fewer vertices where a single virtual particle branches out to two virtual 
particles ("1-2" vertices). This expansion in e is roughly equivalent to the familiar expansion 
in the continuum theory or in quantum field theory in terms of the number of loops. Using 
this expansion, we can numerically evaluate the successive derivatives of the correction term 
at / = 1, and compare these results to experiment. The corrections thus calculated agree 
exactly with our experimental results to within the statistical accuracy of the experimental 
results. This calculation has the additional feature that it is possible to prove that the sum 
over all diagrams which contribute to a given order in e converges. This convergence follows 
from the fact that the sum over all diagrams which contain a fixed number of 1-2 vertices 
is convergent; we prove this for the case of a single 1-2 vertex using an argument which can 
be generalized in a straightforward fashion. In Subsection [7.4| we consider the ring approxi- 
mation from an analytic perspective, and compare this approximation to the closely related 
2-particle BBGKY approximation. Finally, in Subsection |7|.5 we briefly compare the results 
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of the varied approximation methods used in this section. 



7.1 The Short-r Approximation 

Consider the corrections to the eigenvalue A 2 in the Chapman-Enskog analysis of the 1D3P 
lattice gas arising from diagrams of fixed length r. The first few such corrections are easy 
to evaluate by hand. For r = 3, the only diagrams contributing to Eq. Q6.4Q are the two 
diagrams 7\ and T2 shown in Fig. [| These two diagrams shift the eigenvalue A 2 by 

5A 2 = A 2 ( 3 )-A 2 = 9p 2 f(l-f)[(fp) 2 -fp(l-2fp)] 

= - 9 p 3 / 2 (l-/)(l-3M (7.10) 

where we denote by A 2 *- 3 - 1 the eigenvalue of qf in the r = 3 matrix 1 ■. For r = 4, there are 

22 diagrams which contribute to Eq. ( |6.4|) . The 11 diagrams contributing to Y ± are shown 

in Fig. [|; the diagrams contributing to Y + ~ can be generated from these by simply changing 
the directions of the final pair of virtual particles. The correction to A 2 for r = 4 is given by 

5A 2 = A 2 ^ - A 2 ^ = 9p 4 / 2 (l - /)(1 - 3/p) (4 - 3/ - 3f 2 p) . (7.11) 



The complete correction in the r = 4 short-r approximation is given by summing the shifts in 
Eqs. ( 7. 10|) and ( 7.11 ). As r increases, the number of diagrams contributing to 5 A 2 increases 



exponentially, and it rapidly becomes impractical to compute the exact correction, even using 
numerical computing techniques, without some means of simplifying or approximating the 
calculation. We have calculated the corrections to A 2 including diagrams up to r = 5. The 
resulting short-r approximations A 2(r ) are plotted against / and compared to experiment for 
p = 1/4 and 1/2 in Figs. ||] and[l(| respectively. Although the results of this calculation are a 
great improvement over the Boltzmann approximation, and clearly appear to be converging 
to the experimental values as r increases, the oscillatory nature of these approximations (in 
r) is an undesirable feature which makes it difficult to use partial results to put bounds 
on the actual diffusivity. Nonetheless, it is clear from the graphs in Figs. |9] and [U] that 
this formalism includes empirically measurable effects which are completely dropped in the 
Boltzmann approximation. 



7.2 BBGKY Approximations 

We now consider the partial BBGKY approximations for the 1D3P lattice gas. As described 
in Section |^, the fc-particle BBGKY approximation is given by summing over all diagrams 
which have at most k simultaneously correlated virtual particles. For each k, the approxima- 
tion thus consists of an infinite number of diagrams. This approximation can be reduced to 
a finite sum by also limiting the lengths of the diagrams to some maximum size r as in the 
short-r approximation. Because for fixed k, the computational complexity of the summation 
of graphs of length < r grows polynomially in r rather than exponentially as in the short-r 
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Figure 9: Short-r Approximations for p — 1/4 



approximation, it is easier to compute the limit of the set of fc-particle diagrams as r — > oo 
than the complete set of diagrams in this limit. For k = 2, it is possible to prove that this 
infinite sum of diagrams must in fact converge; we derive this result later in this subsection. 
For k > 2, we do not have a complete proof of convergence; however, numerical evidence 
indicates that for each k, the infinite /c-particle BBGKY sum of diagrams is convergent. By 
using methods like those used in the following subsection, it may be possible to prove that 
for each k the BBGKY approximation converges. 

We have used a computer to numerically calculate the limit of the full /c-particle BBGKY 
approximation for certain values of / and p. The algorithm we used was to sum all diagrams 
of length < r on a lattice of width /, then to take the limits as r, I — > oo. As an example, for 
a characteristic pair of values (p, f) = (.25, .5), we have graphed in Fig. |TT] the corrections 
due to 2-particle BBGKY diagrams of length < r on a lattice of size / for all r < 30 and 
for lattice sizes 2, 3, 4, 8, and 16. The curve for I = 8 is indistinguishable from, and hidden 
by, the curve for I = 16. Note that all diagrams of length r are correctly summed as long 
as I > t for all k; for k = 2, however, many diagrams cancel so that larger values of r give 
exact results, such as for I — 2, r = 4. It is clear from this graph not only that the sum 
over diagrams converges rapidly, but also that the major part of the sum arises from the 
contributions of diagrams which are of limited width. We have numerically approximated 
the limits of the fc-particle BBGKY sums for k < 5. The results of this calculation are 
graphed in Figs. [12] and [13|, for the same ranges of values for p, f which were used for the 
short-r approximations in Figs. ^ and 0. 
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Figure 11: Two-Particle BBGKY Corrections for (p, /) = (0.25,0.5) 
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Figure 13: Partial BBGKY Approximations for p = 1/2 

We conclude this subsection with a proof that the 2-particle BBGKY approximation 
is convergent. Although this proof is carried out for the particular 1D3P lattice gas, it 
generalizes easily to any other system. We begin by considering a matrix M on the space 
B2 of 2-particle subsets of B. The matrix element is defined by summing over all 2- 

particle diagrams which begin with outgoing VP's j at an initial vertex and conclude with 
incoming VP's i at a final vertex, and which have the property that all internal vertices have 
a single incoming and outgoing VP. In particular, this means that not only are we restricting 
to 2-particle diagrams, but we are also selecting that subset of diagrams which have no 
interaction vertices other than the initial and final points. Because the internal vertices V l j 
describing propagation of a single correlated quantity are symmetric in i and j, and satisfy 
the conditions 

< V) < 1 

and 

En = E^ = i, (7.12) 

i i 

it is fairly straightforward to see that the matrix elements of M must satisfy the same 
constraints. Intuitively, the above conditions on V 1 - can be interpreted as arising from a 
description of the propagation of a single correlated quantity as a discrete random walk, 
where the direction of the walk at a given time is dependent on the direction of the walk at 
the previous time step, according to the rule that the (virtual) particle will change direction 
with probability g — fp and will continue moving in the same direction with probability 
1 — 2(7. The possible directions of the walk correspond to the lattice vectors eo, e±. With this 
interpretation, the matrix element M~ gives the probability that a pair of virtual particles 
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which begin in state j will collide for the first time in state i. Because for g > the 
probability is 1 that two particles beginning at the same lattice site will eventually collide, 
it follows that condition ( 7.12 ) must hold for the matrix M. The remaining conditions on 



M of symmetry and positivity follow immediately from the random walk interpretation. 

It is now possible to give a complete description of the 2-particle BBGKY approximation 
in terms of the matrix M. The 2-particle BBGKY correction to the transport coefficients of 
any lattice gas can be described in terms of the correction to the J-matrix 

0J j — V lm V j 1V1 npi 

where M is a matrix in the space B 2 defined by summing over all 2-particle diagrams. 
This matrix is defined in a similar fashion to M; however, now we are including diagrams 
which have an arbitrary number of 2-particle interaction vertices V~. By considering these 
2-particle interactions as another matrix V in the space B 2 , we can write a matrix equation 
for M in terms of M and V, 

M = M + MVM. 
Formally, this equation has the solution 

M = M(l-VM)-\ (7.13) 

The matrix V restricted to the space B 2 is easily seen to be another positive definite symmet- 
ric matrix which satisfies Eq. ( |7.12| ). It follows that the matrix VM has the same properties. 



Unfortunately, this appears to lead to a difficulty; namely, it is a consequence of Eq. ( [7.12[) 



that all these matrices have eigenvalues of 1, with associated eigenvectors (1,1,1). This 
means that the matrix 1 — VM has a eigenvalue, and thus has no inverse. However, 
we have an additional factor working in our favor. That is, we are not interested in com- 
puting the complete matrix M. Rather, we are interested in computing the part of that 
matrix which contributes to the renormalization of the eigenvalue A 2 . The correction to this 
eigenvalue due to the 2-particle BBGKY approximation is given by 

5X 2 = 9p 2 f(l - f) [M\ - M\] . (7.14) 

The quantity in brackets can be rewritten in matrix form as 

q 2 Mq 2 = q 2 (M + MVM + MVMVM + ■ • ■) q 2 , 

where in the basis (+,0,—), q 2 = (—1,0,-1-1) is the usual left eigenvector of the J-matrix, 
and q 2 is the corresponding right eigenvector. By the same symmetry arguments which we 
used in Subsection |5j] to prove that the vectors q l are eigenvectors of the renormalized J- 
matrix, it follows that q 2 must be an eigenvector of M with some eigenvalue I. Similarly, q 2 
is an eigenvector of V (considered as a matrix in the space B 2 ). It is straightforward to verify 
from the vertex rules that the eigenvalue of q 2 in this matrix is 1 — 3p + 3pf. Thus, due to 
the fact that q 2 is orthogonal to the eigenvector (1, 1, 1) we avoid the divergence associated 
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with the unit eigenvalue of this vector. In terms of the eigenvalue I it is possible to rewrite 
Eq. ([7.14 ) in the form 

A\2 p 2 /(l ~ f)l 

6X =9 1 Tn q 10 fv 7 - 15 

1-1(1- 3p + 3pf) 

which is manifestly a finite renormalization. We have thus shown that for the 1D3P lattice 
gas, the 2-particle BBGKY correction is finite for any p,f. To check this result, we have 
numerically estimated the matrix M for certain values of p, f and verified that Eq. ( |7. 15|) 



gives a correction which agrees with the numerical results obtained from the general BBGKY 
computer code described above. As an example, for (p, /) = (.25, .5) we get a matrix M 
which is approximately 



M 



( 0.305894 0.293698 0.400407 
0.297612 0.404777 0.297612 
V 0.400407 0.293698 0.305894 



The eigenvalue of q 2 for this matrix is I = —0.094513. The estimated 2-particle BBGKY 
correction to the eigenvalue A 2 is thus 

9/ 

SX 2 = — 0.012550, 

64 - 40/ 

in exact agreement with the results computed numerically and graphed in Fig. [II]. 

A particularly simple example of this formalism arises in the case g = fp = 1/3. In 
this case, the random walk which M is described in terms of is a true random walk, with 
the probability at each time step of each of the 3 possible directions being exactly 1/3, 
independent of the direction of the previous step. It follows immediately from a consideration 
of the matrix M that the eigenvalue I of q 2 is in this case 0, which implies that the shift 
to the diffusivity arising from the 2-particle BBGKY approximation is whenever g — 1/3. 
Note, however, that the higher k BBGKY approximations do not generally vanish in this 



case. For example, in Fig. [13] the 2-particle BBGKY approximation vanishes at the point 
(P) /) = (1/2,2/3); however, the higher k approximations do not vanish and are closer to 
the experimental results. 



7.3 Expansion Around / = 1 

In this subsection we consider the renormalization of A 2 when the particle density / ap- 
proaches 1. In this regime, it is possible to expand the eigenvalue A 2 , and thus the diffusivity, 
in the quantity e = 1 — f . By rewriting the correlation vertex coefficients in terms of e, we 
can ascertain which diagrams contribute to each order in e. The CVC's are shown in terms 
of e and g = fp in Fig. [TJ1 Note that both vertices which take a single incoming virtual par- 
ticle to two virtual particles (1-2 vertices) are proportional to e. Since these vertices are the 
only nonzero vertices which increase the number of virtual particles in a diagram, it follows 
that the set of diagrams which contribute to order e n must be a subset of the diagrams in 
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Figure 14: Correlation Vertex Coefficients in Terms of g, e 

the (n + l)-particle BBGKY approximation. In fact, aside from the interactions described 
by the vertices with 2 and 3 virtual particles in both the incoming and outgoing states, the 
ordering of diagrams in e is equivalent to the loop ordering of diagrams which is commonly 
used in continuum kinetic theory and quantum field theory. Because the sets of diagrams 
which contribute to the corrections for low orders in e are fairly simple (but infinite), we can 
numerically evaluate these partial diagrammatic sums to get a prediction for the low-order 
derivatives of D around / = 1. Because the factor g appears in the 1-1 vertices VI-, it is 
convenient to fix this quantity while evaluating the derivatives of D with respect to e. The 
numerical calculation of the low-order coefficients as a sum over diagrams converges quite 
rapidly. In fact, it can be shown that for each n, the coefficient of e n gives a convergent sum; 
we outline a proof of this fact at the end of this subsection. 

As the simplest example of this type of calculation, consider the set of diagrams which 
contribute to order e in A 2 . Because the correction factor in Eq. (|6.5| ) is itself proportional to 
e, the diagrams which contribute linearly in e cannot contain any internal 1-2 vertices; the 
contribution to the eigenvalue from such diagrams is equivalent to the 2-particle BBGKY 
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Figure 15: Quadratic and Cubic Approximations to A 2 Around / = 1 for g = fp = 0.48 

contribution using the limits of the 2-2 vertices as e — > 0. Using the results of the previous 
subsection, this contribution is finite, and can be computed just as the 2-particle BBGKY 
correction is calculated above. As an example, consider again the case where g = 1/3. In 
this case, we expect the derivative of A 2 at / = 1 to be 0; this result seems to be in agreement 
with experiment. 

We have numerically calculated the first 3 derivatives of A 2 at the point / = 1 for several 
values of the parameter g. The results are completely in agreement with experimental results 
when 1-/<1. As an example, we have graphed the quadratic and cubic approximations 
to 5X 2 = A 2 — A 2 for g = 0.48 in Fig. 15, and compared to experimental data at the point 
/ = 0.96. For comparison, the curves describing the r = 3 and r = 4 short-r approximations 
are also graphed in this region. Note that the region of the graph with / < 0.96 corresponds 
to p > 0.5 and is unphysical. 

We conclude this subsection with an outline of the proof that the sum over diagrams 
contributing to the e n term in A 2 is convergent for any fixed value of n. This proof is 
similar in nature to the proof of convergence for the 2-particle BBGKY approximation in 
the previous subsection, but is slightly more subtle. As mentioned above, to prove the 
desired result it will suffice to show that the contribution from all diagrams with n vertices 
of the 1-2 type gives a convergent sum for all values of n. As the simplest example beyond 
the 2-particle BBGKY approximation, we consider the contribution to the order e 2 term 
arising from diagrams with a single internal 1-2 vertex. Using an analogous notation to that 
used in the previous subsection for the analysis of the 2-particle BBGKY approximation, 
the contribution from the diagrams with a single internal 1-2 vertex can be written as 



5\ 2 = 9p 2 /(l - /) (M 



M~ 



79 



V 



V 




j w v y x % 

(1-MV) _1 M M(l-VM)- 1 Af(l-VM) -1 

Figure 16: Diagrammatic Sketch 



where M is given by 



M 1 - 



[(1-MV)-' 1 ' 
x[(l - VM) 



k M k v\My z 



(7.16) 



In this equation, the variables x and w are summed over all distinct elements of the set of 
pairs of bits at distinct lattice points, modded out by equivalence under translation. The 
elements of this set B\ \ are in 1-1 correspondence with triplets (n,i,j), with n > giving 
the distance between the lattice points of the two bits, and i,j denoting the elements of 
the set B corresponding to the two bits. The variables y, z, and v are similarly summed 
over all elements of the set B\i of triplets (n,i,j) representing a single bit i at a distance 
n from a pair of bits j £ Bj. This set contains elements with n both positive and negative, 
since the lattice point with a single VP can be on either side of the lattice point with two 
VP's; there is a single element of the set -81,2 with n = corresponding to the situation 
where all 3 particles are at the same lattice site. With this constraint on the summation 
variables, the matrices M and V are defined in an analogous fashion to the 2-particle BBGKY 
case. For example, the matrix M y z is the matrix on the space B\ 2 whose entries correspond 
to the probability that a set of 3 particles beginning in state z will first have a collision 
when they are in the state y. Similarly, the matrix M k gives the sum over all 2-particle 
diagrams with no (internal) collisions which begin in the state described by x and end with 
the collision described by the state k. The matrix V x y always describes the product of 
correlation vertex coefficients involved in a collision from a state y to a state x\ when both y 
and x are in B\^, we insist that the two particles remain at the same lattice site. Note that 
if a collision involves virtual particles at more than one lattice point, the distances between 
the vertices must be commensurate in both states for this matrix element to be nonzero. In 
Eq. ( |7.16| ) , the matrix inverses of the form (1 — x)^ 1 should be taken to be shorthand for 
the formal expansion J2i x *, as these matrices generally have a unit eigenvalue. Just as in 
the BBGKY case, it is straightforward to verify that all matrices considered here which are 
square (have both indices taking values in the same space) are symmetric, nonnegative, and 
satisfy condition ( |7.12j) . Eq. ( |7.16| ) can be sketched diagrammatically as in Fig. |R. 

By noting that M and its expression ( |7.16|) are symmetric, an identical argument to that 
in the previous subsection tells us that the contribution to the eigenvalue A 2 can be described 
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by the matrix formula 



SX 2 o p'fi 1 - fr 

[l-/(l-3p + 3p/)]s 



X 



q l [M\V x y My z [{\ - VM)- l ] z v V\M^ q™. (7.17) 



We have thus reduced the problem of proving that Eq. ( |7.16| ) is convergent to the problem 
of proving convergence for Eq. ( |7.17p . There are several key arguments necessary to proving 
the convergence of this remaining sum. The first step is to prove that any matrix element 
of the form 

M* w = V x y My z [(l-VM)- l ] z v V v w 

is convergent, and has an absolute value bounded above by some number S. The second 
step is to argue that the infinite sum 

E M ^ (7-18) 

uieSi,i 

converges for all m and thus can be bounded above by another number A. Once these two 
facts are shown, it follows immediately that total eigenvalue shift is bounded above by 

|,A 1 K 9 [i-;(i-3 P + 3 P /)r A ' 

We will now proceed to prove these two necessary convergence results. We first show 
that there is an upper bound S on the matrix elements M x w . For a fixed value of x, the 
matrix elements V x form a vector in i? 12 which we will refer to as t y . An examination of 
the collision rules ( |6 . 1|) tells us that t has only 18 nonzero components; 8 of magnitude —g 2 £, 
4 of magnitude 2g 2 e, 4 of magnitude —g(l — 2g)e and 2 of magnitude 2g(l — 2g)e. We will 
now consider the result s of multiplying this vector by the matrix M 

s z = t y My z . 

Because the sum of the components of t is 0, the same must be true of the vector s. Fur- 
thermore, if we define a norm on t by 



1*1 = El* 



V I J 



then we can proceed to show that 

\s\<\t\{l-g\l-2gf). (7.19) 

This result follows because there is always at least one set of diagrams of order g 6 (l — 2 g) 2 
which cancel between the positive and negative elements of t (in practice, much cancellation 
occurs; however, we are interested here only in the convergence.) An example of three 
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Figure 17: Cancelling Diagrams 



diagrams giving such a cancellation is given in Fig. [L7|. The vertex factors from these 
diagrams are identical; however, the diagrams connect to components of t y with opposite 
sign. We do not depict the motion of the extra particle, which can be assumed to be constant 
for the three diagrams. From Eq. (|7.19| ), it follows that every matrix element of M is bounded 
above by 

2Apge 



\M X J < E 



g 6 (i-2 9 y 



Thus, we have an upper bound of the desired form for the matrix elements of M. 

Finally, it remains to demonstrate that the infinite sum (|7.18|) is bounded above for all 



m. This sum, however, is simply equal to the expected value of the number of time steps 
necessary for a pair of random walks beginning in state m to collide. To see this, observe 
that every diagram with weight W(T) containing two noncolliding random walks which con- 
tributes to M 1 ^ will contribute to the sum (|7.18|) a total of tW(T) where r is the length of 
the diagram; the factor r appears because the diagram can be chopped in half at any point x 
and will contribute a factor of W(T) for each such division. Given this interpretation, how- 
ever, it is clear that the sum is convergent due to the standard result that in one dimension 
a random walk will return to any point on the lattice in a finite expected time (this result 
is usually stated for random walks without memory; however, a generalization to random 
walks with memory and nonzero bounce probability is straightforward.) In conclusion, we 
have proven both convergence bounds which were needed to demonstrate conclusively the 
convergence of the sum over diagrams which contain a single 1-2 vertex. This implies the 
convergence of the order e 2 term in the diffusivity around / = 1. It is fairly straightforward 
to generalize these arguments to the coefficient of an arbitrary order e n by showing that the 
sum over diagrams with any fixed number of 1-2 vertices is convergent; however, the details 
become correspondingly more complex and are left as an exercise to the reader. Note also, 
that we do not have any reason other than empirical results to believe that the expansion 
of A 2 in e has a nonzero radius of convergence, even though the coefficients themselves are 
proven to be finite. 



7.4 Ring Approximation 

We now consider the ring approximation for the 1D3P lattice gas. In this subsection we study 
the ring approximation from two perspectives. First, we show that the ring approximation 
can be formulated in terms of the M matrix in the same form as the 2-particle BBGKY 
approximation, where the collision matrix V is replaced by a new effective collision matrix. 
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This formalism gives an analytic relationship between the corrections to A 2 from the ring 
approximation and the 2-particle BBGKY approximation. Second, we perform an explicit 
analysis of the ring approximation in Fourier space, and describe the complete contribution 
to the eigenvalue renormalization from this approximation. 

The ring approximation is taken by summing over all independent paths for two separate 
virtual particles to propagate from one point to another, using for each virtual particle the 1- 
1 CVC's as weights on the vertices of the independent paths. The set of diagrams associated 
with this approximation is closely related to the 2-particle BBGKY set of diagrams; however, 
there are two important differences. The first essential difference is that because the two 
VP's are moving independently, there is no constraint dictating that the two VP's cannot 
move along the same lattice vector at some time step r. Thus, diagrams such as Fig 



must be included in this approximation. The second essential difference is that even when 
the two VP's enter a vertex from different directions described by the state i and leave in 
different directions j, the amplitude for such a transition is no longer given by V but rather 
by the collision matrix U with elements 

Ul=5 j i [(l-2g)(l-3g)]+g(l-g). 

Despite these differences, it is still possible to formulate an expression for the ring approxi- 
mation which is identical in form to Eq. ( |7. 14j ) , where the matrix M is no longer defined by 
Eq. (17.131) but rather by 

M = Mil-UMy 1 , 

where the matrix U has elements W~ giving the total amplitude for all processes where a 

pair of particles % come together at some vertex, perhaps travel together for several steps, 
and then separate in directions j at the same or a later vertex. Algebraically, U is given by 

U = U + S-T, 

where the matrix T% gives the amplitude for an incoming pair of particles i to both move in 
direction k according to the separate V 1-1 vertices, and the matrix S 3 gives the total am- 

k 

plitude for a pair of particles both moving in direction k to eventually separate in directions 
j. The matrix T can easily be computed from V, and has matrix elements 

T k 7 =6 k M3g-l)]+g(l-2g). 

The matrix S is slightly more complicated. By symmetry, the elements of this matrix are 
given by 

S 7 k = S j k [s-d]+d, 

for some values of the functions s(g), d(g). From the definitions of these functions, one finds 
that they must satisfy the recursion relations 

s = 2g 2 + (l~2g) 2 s + 2g 2 d, (7.20) 
d = 2g(l-2g) + (l-2g) 2 d + g 2 s + g 2 d. (7.21) 
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Figure 18: Unphysical diagram present in the ring approximation 



These equations have the solution 

-3$ 



30-4' 



(7.22) 



d - ^ < 7 - 23 > 

Plugging these values for s and (i into 5, and computing the eigenvalue w of q 2 with respect 
to the resulting matrix U, we find 

and thus 

^ = 1-90^|. (7.24) 
3#-4 

As a result of this analysis, we have an expression similar to Eq. (|7.15| ) for the eigenvalue 
shift, 

SX* = / /(1 7^ , (7.25) 
1 — Iw 

where w is given by Eq. (|7.24|) . We thus have not only managed to prove that the ring 
approximation for the 1D3P lattice gas is convergent, but we can also derive from Eqs. ( pM5|) 
and ( [7.25| ) an explicit algebraic relationship between the eigenvalue renormalization from the 
ring approximation and that from the 2-particle BBGKY approximation. If we write the 
shift in A 2 due to the ring approximation as 5, and the shift from the 2-particle BBGKY 
approximation as we have the equation 

5(4 -3pf) 

5(4 - 21pf + Tiff) + 9p 2 /(l - /)(4 - 3pf) 

( 5/(l_3 J9 + 3p/ ) + 9p2 /(1 _ /) - ( 7 - 26 ) 

As a specific example, we can plug in the explicit value of 5' = —0.012550 for (p, f) = (.25, .5), 
to find that 

5 = -0.012696. 
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This result is in exact correspondence with numerical calculations. 

To derive the correction due to the ring approximation from first principles, note that 
the weight of any ring diagram is the product of contributions from each of the two virtual 
particles. Denote by y l k (x.,i) the factor contributed by a virtual particle beginning at the 
origin at r = in direction k, and terminating at x at r = t in direction %. In terms of these 
factors, the ring approximation can be written 

oo 

= EE(y*( x .'rt x . t )+»i( x .W x . 4 ))- ( 7 - 27 ) 

Furthermore, because the ring approximation considers only one-point to one-point vertices, 
the y l k (x,t) factors obey the linear, homogeneous dynamical equation, 

y i k (x + c\t + At) = V J y\(x,t). 

Because this equation is linear and homogeneous, it can be solved exactly in Fourier space. 
The result can be expressed as a Fourier series in x and t. When these results for |/ l a (x, t) and 
y J b {^, t), with Fourier summation variables k and k 7 respectively, are inserted into Eq. ( |7.27| ), 
the summation over time is nothing more than a geometric series, and one of the Fourier 
summations yields a Kronecker delta in k and k'. The result for b y is then reduced to 

one summation over k. Passing to the hydrodynamic limit, this becomes a quadrature. For 
a .D-dimensional lattice gas, one obtains at most a D-dimensional quadrature. 

Note that this program can be carried out to compute the ring approximation for any lat- 
tice gas. Since the relationship (|7.26| ) between the ring approximation and the k = 2 BBGKY 
truncation is also generalizable to any lattice gas, it follows that the latter approximation is 
also reduced to quadrature. 

7.5 Comparison of Approximations 

We conclude this section with a brief discussion of the relative effectiveness of the various 
methods used to compute partial diagrammatic sums. We have calculated explicitly in this 
section the corrections to the diffusivity in the 1D3P lattice gas which arise from various 
methods of truncating the complete diagrammatic summation. The short-r approximations 
seem to converge, however require exponential time to compute, so that achieving a high 
degree of accuracy with this method is difficult. The fc-particle BBGKY approximations can 
be calculated in polynomial time, and seem to converge rapidly for each value of A;. It seems 
likely that for most other lattice gases of interest, these approximations will behave similarly, 
and thus that in general the BBGKY approximations will be the more efficient of these two 
methods to achieve a maximal degree of accuracy with a minimum of computation. 

We have also discussed the ring approximation, and shown that the result of this ap- 
proximation is closely related to the 2-particle BBGKY approximation through an algebraic 
relation. Since the ring approximation can be performed analytically for most lattice gases, 
at least in an asymptotic sense, this approximation is generally useful for indicating the 
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convergence properties of the diagrammatic summation for any given lattice gas. In general, 
for lattice gases in 2 dimensions which preserve momentum as well as particle number, the 



ring approximation diverges logarithmically in the lattice size ||35|| . This divergence can be 



reproduced from the diagrammatic formalism directly; this issue will be discussed in more 
detail in a future paper p6 |. 

For the 1D3P lattice gas we have studied in this section, we found that in the vicinity 
of the parameter value / = 1, it is possible to expand the sum over diagrams in a power 
series in 1 — /, giving a perturbation series analogous to the loop expansion in field theory 
or continuum kinetic theory. We showed that each term in this expansion corresponds to a 
convergent sum. In other lattice gases, similar expansions may be useful in calculating the 
effect of renormalization in the vicinity of certain parameter values. When expansions of 
this type are possible, they are generally more accurate than any of the other approximation 
methods. 



Finally, in Fig. [L9| we compare numerically the results of the different approximation 
methods described above to the Boltzmann and experimental calculations of the transport 
coefficients. 
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Figure 19: Comparison of approximation methods for 1D3P 
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8 Conclusions 



In this paper, we have presented a complete kinetic theory of lattice gases, applied it to four 
model lattice gases, and compared the predictions of the theory to experiment for one of 
these models. 

The approach presented in this paper opens up a wide range of possible new work on 
discrete kinetic theory. By applying these techniques to compute deviations from the Boltz- 
mann predictions for commonly used lattice gases, the results of simulations can be more 
accurately interpreted. Lattice gases are currently being used, both in industrial and aca- 
demic settings, for computational fluid dynamics calculations p8| ; to ensure the accuracy 



of these calculations, it is essential to account for the renormalization effects that we have 
studied here. 

In addition to quantitative refinement of lattice gas calculations, the theory presented here 
provides a tool with which to investigate fundamental physical phenomena in nonequilibrium 
statistical systems. In recent years, for example, lattice gases have been used to model many 
different hydrodynamic systems, including reaction-diffusion equations, and other systems 
capable of spontaneous self-organization. It is known || that the Boltzmann approximation 
does not yield accurate results for the transport coefficients of such systems, unless the 
reactants are allowed to diffuse for several steps between reactions in order to artificially 
supress the correlations that develop f39fl . Thus, it is likely that the methods described in 
this paper will be directly applicable to these systems, providing an important correction to 
their theory. More interestingly, these methods will also provide insight into the extremely 
subtle flow and agglomeration of interparticle correlations in pattern-forming lattice gases, 
and hence into the dynamical basis of self-organization. 
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A Glossary of Notation 



In this appendix, we list all the important symbols used in this paper, giving the page number 
where they were first used (if appropriate) and a brief description (if appropriate). 
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DESCRIPTION 

Used in place of an index to indicate possible func- 
tional dependence on all values of an index 
Empty set 

Placed over a set to indicate that it contains 

strictly more than one element 

Equivalence relation between two states 

Used to indicate subsets of the set B of bits in the 

entire lattice 

Parameters of Fermi-Dirac equilibrium 

Set of bit numbers for all particles in /3 at site x 

Connected correlation function 

Fermi connection 

Generalized Fermi connection 

Kronecker delta of two bits 

Perturbation expansion parameter 

1 — / for the 1D3P lattice gas 

Eigenvalue of J 4 - 

Notation for {i(a\ ),..., i{o>j)} 

Used to indicate subsets of the set B of bits at a 

single lattice site 

Set of outgoing virtual particles in diagram T at 
timestep r 

Set of incoming virtual particles in diagram T at 
timestep r 

Set of all partitions of set a 

Number of outgoing virtual particles in diagram 

T at timestep r 

Used to enumerate all the bits on the lattice 

Bit number on the lattice of the ith bit at site x 

Microscopic transition matrix element 

Ensemble-averaged transition matrix element 

Advection coefficients in hydrodynamic equation 

Advection operator 

The set of n bits at a site 

The set of iV bits on the lattice 

Characteristic lattice spacing 

Microscopic collision operator 

Lattice vector along which particles can move 
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Ensemble averaged collision operator 

£th order, ensemble-averaged collision operator 

Spatial dimension of lattice 

Diffusion coefficients in hydrodynamic equation 

Dimensionless lattice vectors 

Functional dependence of multipoint means on 

connected correlation functions 

Functional dependence of connected correlation 

functions on multipoint means 

Fermi metric 

Generalized Fermi Metric 

Set of hydrodynamic modes, {1, . . . , n c } 

Generalized Kronecker delta 

Used to enumerate all the bits at a given site 

The bit number (at its site) of the ath bit on the 

lattice 

Jacobian matrix of collision operator at equilib- 
rium 

Renormalized Jacobian of collision operator 
Exact propagator for one-point means 
Coefficient of dependence of c l on n^w? 
Length of diagram T 
Set of kinetic modes, {n c + 1, . . . , n} 
Collision operator for multipoint means 
Collision operator for connected correlation func- 
tions 
Lattice 

Subset of points in L which contain at least one 

particle in set ft 

Number of bits per site 

ith bit at site x at timestep t 

The value of the ath bit on the lattice at timestep 

t 

Sometimes used to denote stochastic bits, when 
no ambiguity will arise. 
Number of conserved quantities 
Total number of bits on the lattice 
Multipoint mean 
Ensemble average of n*(x, t) 
ith order, ensemble-averaged distribution func- 
tion 

Multipoint probability distribution 

Value of /ith conserved quantity at position x at 

timestep t 
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Source term in hydrodynamic equation 
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Factor contributed by one vertex, at site x, to /C 
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Weight of diagram T 
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Position of a site on lattice L 


x(a) 




The site of the ath bit on the lattice 
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